#include "mcmcRegionModel.h"
#include "mcmcMap.h"
#include "mcEdgeOp.h"

double scale_factor	=	17.5;

#define COMPLEXTY_SINGLEGAU		0.00
#define COMPLEXTY_MIXTUREGAU	0.30
#define COMPLEXTY_2DCOLOR		log(1.0)
#define COMPLEXTY_3DCOLOR		log(1.0)
#define COMPLEXTY_MIXTURE3DCOLOR	log(100.0)
#define COMPLEXTY_TEXTURE		-1.5
#define COMPLEXTY_INTENSITY		log(1.0)
#define COMPLEXTY_MIXTUREINTENSITY log(1.15)
#define COMPLEXTY_TEXTUREINTENSITY	0.72
#define COMPLEXTY_DELTA			0.25
// parameters for big zebra image
//#define COMPLEXTY_TEXTUREINTENSITY	1.1
//#define COMPLEXTY_DELTA			0.4

#define COMPLEXTY_SURFACE		0.07
#define COMPLEXTY_COLOR_SURFACE	0.25

#define BASE_PIXEL_SHORTEN   2.0;

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
MCMC_Model::MCMC_Model()
{
	prespFilter			= NULL;
}

MCMC_Model::~MCMC_Model()
{
	prespFilter = NULL;
}

MCMC_Model& MCMC_Model::operator=(const MCMC_Model&x)
{
	if (this!=&x)
	{
		vc_vote			= x.vc_vote;
		crying_voice	= x.crying_voice;
		log_pos_prob	= x.log_pos_prob;
		nsize			= x.nsize;
		log_lik			= x.log_lik;
		prespFilter		= x.prespFilter;
	}
	return *this;
}

void MCMC_Model::SaveToFile(ofstream &strm)
{
	strm.write((char*)&type, sizeof(type));

	vc_vote >> strm;

	strm.write((char*)&crying_voice, sizeof(crying_voice));
	strm.write((char*)&log_pos_prob, sizeof(log_pos_prob));
	strm.write((char*)&nsize, sizeof(nsize));
	strm.write((char*)&log_lik, sizeof(log_lik));
}

void MCMC_Model::ReadFromFile(ifstream &strm)
{

	// type is read before this
	vc_vote << strm;

	strm.read((char*)&crying_voice, sizeof(crying_voice));
	strm.read((char*)&log_pos_prob, sizeof(log_pos_prob));
	strm.read((char*)&nsize, sizeof(nsize));
	strm.read((char*)&log_lik, sizeof(log_lik));
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

// constructions & destructions
MCMC_RegionModel::MCMC_RegionModel()
:MCMC_Model()
{
	pmx_label	= NULL;
	pparticle	= NULL;
	pmx_taken	= NULL;

	plist_domain = NULL;
}

MCMC_RegionModel::~MCMC_RegionModel()
{
	Free();		// free all the allocated space
}

// free allocated space
void MCMC_RegionModel::Free(void)
{
	pmx_label   = NULL;
	pmx_taken	= NULL;
	pparticle	= NULL;

	plist_domain = NULL;
}

// operators
MCMC_RegionModel& MCMC_RegionModel::operator=(const MCMC_RegionModel &x)
{
	if (this!=&x)
	{
		*((MCMC_Model*) this) = *((MCMC_Model*) &x);

		pparticle	= x.pparticle;
		pmx_label	= x.pmx_label;
		pmx_taken	= x.pmx_taken;
		bound		= x.bound;
	}
	return *this;
}

void MCMC_RegionModel::SaveToFile(ofstream &strm)
{
	MCMC_Model::SaveToFile(strm);
}

void MCMC_RegionModel::ReadFromFile(ifstream &strm)
{
	MCMC_Model::ReadFromFile(strm);
}

// operations
bool MCMC_RegionModel::Set(mcmcMatrix<int> &mask_label, irmm_Rect &bound1, MCMC_Particle &particle, mcmcMatrix<int> *pmx_taken1)
{
	// free allocated space
	Free();
	// assign values for the preparation of mean and variance computations
	pmx_label	= &mask_label;
	pmx_taken	= pmx_taken1;
	pparticle	= &particle;
	bound = bound1;

	if (pparticle != NULL)
		vc_vote.SetRowVector(pparticle->Num());

	return true;
}

bool MCMC_RegionModel::Set(LinkList<PMCMC_PixelSet> &list_domain, MCMC_Particle &particle, mcmcMatrix<int> *pmx_taken1)
{
	// free allocated space
	Free();

	plist_domain = &list_domain;

	pmx_taken	= pmx_taken1;
	pparticle	= &particle;

	if (pparticle != NULL)
		vc_vote.SetRowVector(pparticle->Num());

	return true;
}

mcmcVector<float> MCMC_RegionModel::ComputeVote(const int nlabel)
{
	vc_vote.InitValue(0.0f);

	if (plist_domain == NULL)
		for (int j=bound.top; j<=bound.bottom; j++)
			for (int i=bound.left; i<=bound.right; i++)
			{
				if (pmx_label->Data(j,i)>0 && (*pmx_label)(j,i)==nlabel && ConsiderThisPixel(j,i))
					pparticle->AddClusterVote(vc_vote,j,i);
			}
	else
	{
		bool			bgo,bgo_domain;
		int				i,j;
		PMCMC_PixelSet	ppixel_set;

		bgo_domain = plist_domain->GetFirst(ppixel_set);
		while (bgo_domain)
		{
			if (ppixel_set != NULL)
			{
				bgo = ppixel_set->GetFirst(j, i);
				while (bgo)
				{
					if (ConsiderThisPixel(j,i))
						pparticle->AddClusterVote(vc_vote,j,i);

					bgo = ppixel_set->GetNext(j, i);
				}
			}
			bgo_domain = plist_domain->GetNext(ppixel_set);
		}
	}
	return vc_vote;
}

mcmcVector<double> MCMC_RegionModel::DiffuseTo(MCMC_Particle &particle, mcmcVector<double> &vc_mean, const int ipar)
{
	int				iparticle;
	mcmcVector<double>	vc_result;

	// find the particle to which it belongs
	if (ipar<0)
		iparticle = particle.ClosestParticle(vc_mean);
	else
		iparticle = ipar;
	vc_result = vc_mean;

	particle.Diffuse(vc_result,iparticle);

	return vc_result;
}


void MCMC_RegionModel::RobustVariance(mcmcMatrix<double> &mx_var)
{
	double ddeti = mx_var.Det();
	if (ddeti<1e-2)
	{
		mx_var.Identity();
		mx_var = mx_var*pow(10.0,1/3.0);
	}
	else
	{
		if (ddeti<10.0)
		{
			mx_var.Identity();
			ddeti = 1.0;
			mx_var = mx_var*pow(10.0/ddeti,1/((double)(mx_var.rows())));
		}
		ddeti = mx_var.Det();
	}
}

void MCMC_RegionModel::RobustVariance(double &dvar)
{
	if (dvar<2.0)
		dvar = 2.0;
}

double MCMC_RegionModel::RobustLoglik(const double dvalue)
{
	// if this is changed, the corresponding function in "mcmcLayerCurve" needs to be
	// changed as well
	if (dvalue>60.0/2.0)
		return 60.0/2.0;
	else
		return dvalue;
}

void MCMC_RegionModel::GetXYDirections(mcmcVector<double> &vc_X, mcmcVector<double> &vc_Y, mcmcVector<double> &vc_N)
{
	double din,dtemp;

	vc_X.SetRowVector(3);
	vc_Y.SetRowVector(3);

	din = vc_N.Norm(2.0);
	if (din>0.0)
	{
		vc_N = vc_N/din;
		// X=(0, C/sqrt(B2+C2), -B/sqrt(B2+C2))
		// Y=(1-A2, -AB, -AC)
		int i;
		for (i=0; i<3; i++)
			if (fabs(vc_N(i))<0.9)
				break;
		assert(i>=0 && i<3);
		switch (i)
		{
		case 0:
			dtemp = sqrt(pow(vc_N(1),2.0)+pow(vc_N(2),2.0));
			vc_X(0) = 0.0;
			vc_X(1) = vc_N(2)/dtemp;
			vc_X(2) = -vc_N(1)/dtemp;
			vc_Y(0) = 1-pow(vc_N(0),2.0);
			vc_Y(1) = -vc_N(0)*vc_N(1);
			vc_Y(2) = -vc_N(0)*vc_N(2);
			break;

		case 1:
			dtemp = sqrt(pow(vc_N(0),2.0)+pow(vc_N(2),2.0));
			vc_X(0) = vc_N(2)/dtemp;
			vc_X(1) = 0.0;
			vc_X(2) = -vc_N(0)/dtemp;
			vc_Y(0) = -vc_N(1)*vc_N(0);
			vc_Y(1) = 1-pow(vc_N(1),2.0);
			vc_Y(2) = -vc_N(1)*vc_N(2);
			break;

		case 2:
			dtemp = sqrt(pow(vc_N(0),2.0)+pow(vc_N(1),2.0));
			vc_X(0) = vc_N(1)/dtemp;
			vc_X(1) = -vc_N(0)/dtemp;
			vc_X(2) = 0.0;
			vc_Y(0) = -vc_N(2)*vc_N(0);
			vc_Y(1) = -vc_N(2)*vc_N(1);
			vc_Y(2) = 1-pow(vc_N(2),2.0);
			break;
		}
		vc_X = vc_X/vc_X.Norm(2.0);
		vc_Y = vc_Y/vc_Y.Norm(2.0);
	}
}

double MCMC_RegionModel::Partition4Gaussian(const double mean, const double var)
{
	double x,eng,z,c,temp;

	c = sqrt(2.0*PI*var);
	z = 0.0;
	x = 0.0;
	while (x<=INTENSITY_HIGH+INTENSITY_DELTA-1)
	{
		temp = fabs(x-mean);
		eng = temp*temp/(2.0*var);
		z += ProbFromEnergy(eng)/c;

		x = x+1.0;
	}
	return z;
}


bool MCMC_RegionModel::Mean(const int nlabel)
{
	if (plist_domain == NULL)
		return MeanFromLabel(nlabel);
	else
		return MeanFromPixelSet();
}

bool MCMC_RegionModel::Variance(const int nlabel)
{
	if (plist_domain == NULL)
		return VarianceFromLabel(nlabel);
	else
		return VarianceFromPixelSet();
}

bool MCMC_RegionModel::LogLikelihood(const int nlabel)
{
	if (plist_domain == NULL)
		return LogLikelihoodFromLabel(nlabel);
	else
		return LogLikelihoodFromPixelSet();
}

irmm_Rect MCMC_RegionModel::GetBoundFromDomainList(void)
{
	irmm_Rect	bound_temp;
	bool		bgo_domain;
	MCMC_PixelSet	*ppixel_set;

	bound.Set(9999,9999,-9999,-9999);
	if (plist_domain == NULL)
		return bound;

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bound_temp = ppixel_set->Bound();

			bound.SetBound(bound_temp.left, bound_temp.top);
			bound.SetBound(bound_temp.right, bound_temp.bottom);
		}

		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	return bound;
}


///////////////////////////////////////////////////////
//////////////////////////////////////////////////////

// construction & destruction
MCMC_Model_Texture::MCMC_Model_Texture()
:MCMC_RegionModel()
{
	complexity_prior = COMPLEXTY_TEXTURE;

	type			= MODEL_TEXTURE;

	ktexture		= -1;
}

MCMC_Model_Texture::~MCMC_Model_Texture()
{
	Free();
}

void MCMC_Model_Texture::Free(void)
{
	MCMC_RegionModel::Free();

	ktexture = -1;
}

// operators
MCMC_Model_Texture& MCMC_Model_Texture::operator=(const MCMC_Model_Texture &x)
{
	if (this!=&x)
	{
		*((MCMC_RegionModel*) this) = *((MCMC_RegionModel*) &x);

		vc_mean = x.vc_mean;
		ktexture = x.ktexture;
	}
	return *this;
}

void MCMC_Model_Texture::SaveToFile(ofstream &strm)
{
	MCMC_RegionModel::SaveToFile(strm);

	strm.write((char*)&ktexture, sizeof(ktexture));
	vc_mean >> strm;
}

void MCMC_Model_Texture::ReadFromFile(ifstream &strm)
{
	MCMC_RegionModel::ReadFromFile(strm);

	strm.read((char*)&ktexture, sizeof(ktexture));
	vc_mean << strm;
}

// operations
MCMC_Model* MCMC_Model_Texture::CopyTo(void)
{
	MCMC_Model_Texture *preturn;

	preturn = new MCMC_Model_Texture();
	(*preturn) = (*this);

	return ((MCMC_Model*)preturn);
}

bool MCMC_Model_Texture::Mean(const int nlabel)
{
	//int 	i,j;
	int i, j;

	nsize = 0;
	vc_mean.SetRowVector(3);
	vc_mean.InitValue(0.0);


	// compute the mean value
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && (*pmx_label)(j,i)==nlabel && ConsiderThisPixel(j,i))
			{
				nsize++;

				vc_mean(0) += (double)prespFilter->image_original.R(j,i);
				vc_mean(1) += (double)prespFilter->image_original.G(j,i);
				vc_mean(2) += (double)prespFilter->image_original.B(j,i);
			}
		}
	if (nsize > 0)
		vc_mean = vc_mean/((double)nsize);

	if (nsize>0)	return true;
	else			return false;
}

// log likelyhood
bool MCMC_Model_Texture::LogLikelihood(const int nlabel)
{
	int 		i,j;
	mcmcVector<double>	I_u(3);
	double			dis_uc;

	// initialization
	nsize = 0;
	vc_vote.InitValue(0);
	crying_voice = 0.0;
	log_pos_prob = 0.0;

	MCMC_Particle_Texture *pparticle_ple;
	pparticle_ple = (MCMC_Particle_Texture*)pparticle;

	// use a robust statistics method
	int	num_big,num_ok;
	double	mean_ok;

	num_big = 0;
	num_ok	= 0;
	mean_ok	= 0.0;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && (*pmx_label)(j,i)==nlabel && ConsiderThisPixel(j,i))
			{
				if (pparticle_ple->p_mx_energymap[ktexture](j,i)>50.0)
					num_big++;
				else
				{
					log_pos_prob += pparticle_ple->p_mx_energymap[ktexture](j,i);
					mean_ok +=pparticle_ple->p_mx_energymap[ktexture](j,i);
					num_ok++;
				}

				log_pos_prob += pparticle_ple->mx_jacobs(j,i);

				pparticle->AddClusterVote(vc_vote,j,i);

				// compute the crying voice of the regions
				dis_uc = RobustCryingVoice(pparticle_ple->p_mx_energymap[ktexture](j,i)+Complexity());
				crying_voice = crying_voice+dis_uc;

				nsize++;
			}
		}

	mean_ok = mean_ok/((double)num_ok);
	// use a robust statistics method
	if (num_big<150 && ((double)num_big)/((double)nsize)<0.08)
		log_pos_prob += mean_ok*num_big;
	else
		log_pos_prob += 50.0*num_big;

	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	log_lik = log_pos_prob;

	log_pos_prob += pow((double)nsize, COMPLEXITY_POW)*scale_factor+(double)nsize*Complexity();
	crying_voice = pow(crying_voice, 0.6);

	if (nsize>0)	return true;
	else			return false;
}

bool MCMC_Model_Texture::SynthesizeSampling(mcmcVector<double> &vc_a, const int row, const int col)
{
	MCMC_Particle_ColorTexture *pparticle_ple;
	pparticle_ple = (MCMC_Particle_ColorTexture*)pparticle;

	vc_a(0) = pparticle_ple->p_mx_syn[ktexture].R(row,col);
	vc_a(1) = pparticle_ple->p_mx_syn[ktexture].G(row,col);
	vc_a(2) = pparticle_ple->p_mx_syn[ktexture].B(row,col);

	return true;
}

bool MCMC_Model_Texture::Reconstruct(mcmcVector<double> &vc_a, const int row, const int col)
{
	// this is not a legal function call
	return SynthesizeSampling(vc_a,row,col);
}

// get the log of the posterior probability at (row,col)
double MCMC_Model_Texture::GetLogPostProb(const int row, const int col, mcmcMatrix<int> *pmx_occupied)
{
	MCMC_Particle_Texture *pparticle_ple;
	pparticle_ple = (MCMC_Particle_Texture*)pparticle;

	int		i,j;
	double	dis,count,dis_temp;

	dis = 0.0;
	count = 0.0;
	for (j=row-2; j<=row+2; j++)
		for (i=col-2; i<=col+2; i++)
			if (pparticle_ple->p_mx_energymap[ktexture].Valid(j,i))
			{
				dis_temp = pparticle_ple->p_mx_energymap[ktexture](j,i)+pparticle_ple->mx_jacobs(j,i);
				if (dis_temp > 100.0)	dis_temp = 100.0;

				dis += dis_temp;

				count = count+1.0;
			}
	if (count>0)
		dis = dis/count;

	return dis;
}

// for the computation of current model only
void MCMC_Model_Texture::DiffuseComputedMean(const int nlabel)
{
	// compute all the votes
	ComputeVote(nlabel);
	
	// choose a particle
	ktexture = pparticle->ChooseParticle(vc_vote);
}


void MCMC_Model_Texture::OutputInfo(char *info)
{
	if (info!=NULL)
		sprintf(info, "color texture (%d)-%.2f", ktexture,log_pos_prob);
}

double MCMC_Model_Texture::GetPdf(void)
{
	mcmcVector<double> vc_k(1);

	vc_k(0) = (double)ktexture;
	assert(pparticle!=NULL);
	return pparticle->GetPdf(vc_vote,vc_k);
}

void MCMC_Model_Texture::GetDirectValue(mcmcVector<double> &vc, const int row_pixel,const int col_pixel)
{
	vc(0) = prespFilter->image_original.R(row_pixel,col_pixel);
	vc(1) = prespFilter->image_original.G(row_pixel,col_pixel);
	vc(2) = prespFilter->image_original.B(row_pixel,col_pixel);
}

double MCMC_Model_Texture::GetDirectDis(mcmcVector<double> &vc_pixel_rgb)
{
	double dtemp=0.0;

	for (int i=0; i<vc_pixel_rgb.size(); i++)
		dtemp += pow(vc_pixel_rgb(i)-vc_mean(i),2.0);
	return dtemp;
}


////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
//        Models for grey level images


//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
MCMC_Intensity_Gaussian::MCMC_Intensity_Gaussian()
:MCMC_RegionModel()
{
	complexity_prior = COMPLEXTY_INTENSITY;

	type = MODEL_INTENSITY_GAUSS;

}

MCMC_Intensity_Gaussian::~MCMC_Intensity_Gaussian()
{
	Free();
}

void MCMC_Intensity_Gaussian::Free(void)
{
}

// operators
MCMC_Intensity_Gaussian& MCMC_Intensity_Gaussian::operator=(const MCMC_Intensity_Gaussian &x)
{
	if (this!=&x)
	{
		*((MCMC_RegionModel*) this) = *((MCMC_RegionModel*) &x);

		intensity_mean	= x.intensity_mean;
		intensity_var	= x.intensity_var;
		bin_energy		= x.bin_energy;
	}
	return *this;
}

int MCMC_Intensity_Gaussian::TotalCodingBits(void)
{
	int bits;

	bits = sizeof(uchar)+		// mean
		   sizeof(float);		// variance

	return bits;
}

void MCMC_Intensity_Gaussian::SaveToFile(ofstream &strm)
{
	MCMC_RegionModel::SaveToFile(strm);

	strm.write((char*)&intensity_mean, sizeof(intensity_mean));
	strm.write((char*)&intensity_var, sizeof(intensity_var));
	double Z=0.0;
	strm.write((char*)&Z, sizeof(Z));
}

void MCMC_Intensity_Gaussian::ReadFromFile(ifstream &strm)
{
	MCMC_RegionModel::ReadFromFile(strm);

	strm.read((char*)&intensity_mean, sizeof(intensity_mean));
	strm.read((char*)&intensity_var, sizeof(intensity_var));

	double Z;
	strm.read((char*)&Z, sizeof(Z));
}

MCMC_Model* MCMC_Intensity_Gaussian::CopyTo(void)
{
	MCMC_Intensity_Gaussian *preturn;

	preturn = new MCMC_Intensity_Gaussian();
	(*preturn) = (*this);

	return ((MCMC_Model*)preturn);
}

// compute mean and variance for the regions
bool MCMC_Intensity_Gaussian::MeanFromLabel(const int nlabel)
{
	// compute the intensity means
	int 	i,j;

	intensity_mean = 0.0;
	nsize = 0;

	// compute the major direction
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				// get intensity
				intensity_mean += prespFilter->ImageIntensity(j,i);

				nsize++;
			}
		}

	if (nsize > 0)
		intensity_mean = intensity_mean/((double)nsize);

	// robust mean
	double mean_temp=0.0;
	double size_temp=0.0;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				// get intensity
				if (fabs(prespFilter->ImageIntensity(j,i)-intensity_mean)<5.0)
				{
					mean_temp += prespFilter->ImageIntensity(j,i);
					size_temp += 1.0;
				}
			}
		}
	if (nsize>300 && size_temp/nsize>0.90)
		intensity_mean = mean_temp/size_temp;

	if (nsize>0)	return true;
	else			return false;
}

bool MCMC_Intensity_Gaussian::MeanFromPixelSet(void)
{
	// compute the intensity means
	int		i,j;
	bool	bgo,bgo_domain;
	PMCMC_PixelSet ppixel_set;

	intensity_mean	= 0.0;
	nsize			= 0;

	assert(plist_domain != NULL);
	
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// get intensity
					intensity_mean += prespFilter->ImageIntensity(j,i);
					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}
	if (nsize>0)
		intensity_mean = intensity_mean/((double)nsize);

	// robust mean
	double mean_temp=0,size_temp=0;

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// get intensity
					if (fabs(prespFilter->ImageIntensity(j,i)-intensity_mean)<5.0)
					{
						mean_temp += prespFilter->ImageIntensity(j,i);
						size_temp += 1.0;
					}
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (nsize>300 && size_temp/nsize>0.90)
		intensity_mean = mean_temp/size_temp;

	if (nsize>0)	return true;
	else			return false;
}

bool MCMC_Intensity_Gaussian::VarianceFromLabel(const int nlabel)
{
	int 	i,j;
	double		din,num_small,small_var,num_count;

	// initialization
	nsize = 0;
	intensity_var = 0.0;
	num_small = 0.0;
	small_var = 0.0;
	num_count = 0.0;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				// projection
				din = prespFilter->ImageIntensity(j,i);

				// for robust variance
				if (fabs(din-intensity_mean)<5.0)
				{
					num_small += 1.0;
					small_var += pow(din-intensity_mean,2.0);
					num_count += 1.0;
				}

				// variance of intensity
				intensity_var += pow(din-intensity_mean,2.0);
				nsize++;
			}
		}

	if (nsize>1)
	{
		if (nsize>300 && num_small/((double)nsize)>0.90)
			intensity_var = small_var/num_count;
		else
			intensity_var = intensity_var/((double)nsize-1);
	}
	else
		intensity_var = 1.0;

	if (intensity_var<5.0)
		intensity_var = 5.0;

	// compute quantized energy
	Histogram	hist_temp;
	double		value=INTENSITY_LOW,prob,dis;

	bin_energy.Set(INTENSITY_LOW,INTENSITY_HIGH,INTENSITY_DELTA);
	hist_temp.Set(INTENSITY_LOW,INTENSITY_HIGH,INTENSITY_DELTA);

	while (value<INTENSITY_HIGH+INTENSITY_DELTA)
	{
		dis = pow(value-intensity_mean,2.0)/intensity_var/2.0;
		prob = ProbFromEnergy(RobustLoglik(dis)+log(sqrt(2.0*PI))+log(intensity_var)*0.5);
		hist_temp.AddBin(value,prob);
		value += 1.0;
	}

	for (i=0; i<hist_temp.count; i++)
	{
		if (hist_temp.pbins[i]>0.0)
			bin_energy.pbins[i] = -log(hist_temp.pbins[i]);
		else
			bin_energy.pbins[i] = 60;
	}

	return true;
}

bool MCMC_Intensity_Gaussian::VarianceFromPixelSet(void)
{
	int		i,j;
	double	din,num_small,small_var,num_count;
	bool	bgo;

	// initialization
	nsize			= 0;
	intensity_var	= 0.0;
	num_small		= 0.0;
	small_var		= 0.0;
	num_count		= 0.0;

	assert(plist_domain != NULL);

	bool			bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// projection
					din = prespFilter->ImageIntensity(j,i);

					// for robust variance
					if (fabs(din-intensity_mean)<5.0)
					{
						num_small += 1.0;
						small_var += pow(din-intensity_mean, 2.0);
						num_count += 1.0;
					}

					// variance of intensity
					intensity_var += pow(din-intensity_mean, 2.0);

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (nsize > 1)
	{
		if (nsize>300 && num_small/((double)nsize)>0.90)
			intensity_var = small_var/num_count;
		else
			intensity_var = intensity_var/((double)nsize-1);
	}
	else
		intensity_var = 1.0;

	if (intensity_var < 5.0)
		intensity_var = 5.0;
	else if (intensity_var > 1e5)
		intensity_var = 1e5;

	// compute quantized energy
	Histogram	hist_temp;
	double		value=INTENSITY_LOW,prob,dis;

	bin_energy.Set(INTENSITY_LOW,INTENSITY_HIGH,INTENSITY_DELTA);
	hist_temp.Set(INTENSITY_LOW,INTENSITY_HIGH,INTENSITY_DELTA);

	while (value<INTENSITY_HIGH+INTENSITY_DELTA)
	{
		dis = pow(value-intensity_mean,2.0)/intensity_var/2.0;
		prob = ProbFromEnergy(RobustLoglik(dis)+log(sqrt(2.0*PI))+log(intensity_var)*0.5);
		hist_temp.AddBin(value,prob);
		value += 1.0;
	}

	for (i=0; i<hist_temp.count; i++)
	{
		if (hist_temp.pbins[i]>0.0)
			bin_energy.pbins[i] = -log(hist_temp.pbins[i]);
		else
			bin_energy.pbins[i] = 60;
	}
	return true;
}


bool MCMC_Intensity_Gaussian::LogLikelihoodFromLabel(const int nlabel)
{
	int 	i,j;
	double		dis_uc,din;

	// initialization
	nsize	= 0;
	log_pos_prob = 0.0;
	vc_vote.InitValue(0);
	crying_voice = 0.0;

	// the real size of the region
	double	real_size=0.0;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && (*pmx_label)(j,i)==nlabel)
			{
				real_size += 1.0;

				if (ConsiderThisPixel(j,i))
				{
					// LogLikelihood in intensity
					din = prespFilter->ImageRegionLayerIntensity(j,i);

					// if this pixel is covered by a base then the energy should be shortened
					log_pos_prob += bin_energy.GetBin(din);

					pparticle->AddClusterVote(vc_vote,j,i);

					// compute the crying voice of this region
					dis_uc = RobustCryingVoice(bin_energy.GetBin(din)+Complexity());
					crying_voice += dis_uc;

					nsize++;
				}
			}
		}

	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	// the third direction is treated as an independent variable
	log_lik = log_pos_prob;

	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;

	log_pos_prob += ((double)nsize)*Complexity();

	crying_voice = pow(crying_voice,0.6);

	return true;
}

bool MCMC_Intensity_Gaussian::LogLikelihoodFromPixelSet(void)
{
	int		i,j;
	double	dis_uc,din;

	// initialization
	nsize	= 0;
	log_pos_prob = 0.0;
	vc_vote.InitValue(0);
	crying_voice = 0.0;

	// the real size of the region
	double	real_size,dtemp;
	real_size = 0.0;

	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			ppixel_set->energy = 0.0;

			// get each pixel
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				real_size += 1.0;

				if (ConsiderThisPixel(j,i))
				{
					// LogLikelihood in intensity
					din = prespFilter->ImageRegionLayerIntensity(j,i);

					// if this pixel is covered by a base then the energy should be shortened
					dtemp = bin_energy.GetBin(din);
					ppixel_set->energy += dtemp;

					log_pos_prob += dtemp;

					pparticle->AddClusterVote(vc_vote,j,i);
					// compute the crying voice of this region
					dis_uc = RobustCryingVoice(bin_energy.GetBin(din)+Complexity());
					crying_voice += dis_uc;

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	// the third direction is treated as an independent variable
	log_lik = log_pos_prob;
	log_pos_prob += pow((double)nsize, COMPLEXITY_POW)*scale_factor;
	log_pos_prob += ((double)nsize)*Complexity();
	crying_voice = pow(crying_voice, 0.6);

	return true;
}

void MCMC_Intensity_Gaussian::UpdateLogLikelihood(const int nlabel, MCMC_PixelSet &list_pixels_old, MCMC_PixelSet &list_pixels_new)
{
	int 	i,j;
	double		din,dtemp,dcomplexity_temp;
	bool		bgo;

	if (list_pixels_old.Num()<=0 && list_pixels_new.Num()<=0)
		return;

	log_pos_prob -= pow((double)nsize,COMPLEXITY_POW)*scale_factor;

	// delete the old pixels
	dcomplexity_temp = Complexity();
	bgo = list_pixels_old.GetFirst(j, i);
	while (bgo)
	{
		if (pmx_label->Data(j,i) == nlabel)
		{
			// LogLikelihood in intensity
			din = prespFilter->ImageRegionLayerIntensity(j,i);
			dtemp = bin_energy.GetBin(din);

			log_pos_prob += dtemp+dcomplexity_temp;
			log_lik += dtemp;

			nsize++;
		}
		bgo = list_pixels_old.GetNext(j,i);
	}

	bgo = list_pixels_new.GetFirst(j, i);
	while (bgo)
	{
		if (pmx_label->Data(j,i) == nlabel)
		{
			// LogLikelihood in intensity
			din = prespFilter->ImageRegionLayerIntensity(j,i);
			// if this pixel is covered by a base then the energy should be shortened
			dtemp = bin_energy.GetBin(din);

			log_pos_prob = log_pos_prob-(dtemp+dcomplexity_temp);
			log_lik -= dtemp;

			nsize--;
		}
		bgo = list_pixels_new.GetNext(j, i);
	}

	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;
}


void MCMC_Intensity_Gaussian::GetDirectValue(mcmcVector<double> &vc, const int row_pixel,const int col_pixel)
{
	vc(0) = prespFilter->ImageIntensity(row_pixel,col_pixel);
}

double MCMC_Intensity_Gaussian::GetDirectDis(mcmcVector<double> &vc_intensity)
{
	double	din;

	din = vc_intensity(0);

	return (intensity_mean-din)*(intensity_mean-din); // give lower prevelige over intensity
}

void MCMC_Intensity_Gaussian::DiffuseComputedMean(const int nlabel)
{
	mcmcVector<double> vc_temp(1);
	int ipar;

	// compute all the votes
	ComputeVote(nlabel);
	// choose a particle
	ipar = pparticle->ChooseParticle(vc_vote);

	// diffusion
	vc_temp(0) = intensity_mean;
	vc_temp = DiffuseTo(*pparticle,vc_temp,ipar);

	// intensity
	intensity_mean = vc_temp(0);
}

double MCMC_Intensity_Gaussian::GetPdf(void)
{
	mcmcVector<double> vc_temp(1);
	assert(pparticle!=NULL);

	vc_temp(0) = intensity_mean;
	return pparticle->GetPdf(vc_vote,vc_temp);
}

void MCMC_Intensity_Gaussian::PrepreForSynthesize(void)
{
//	rf_i.Set(intensity_mean,sqrt(min(intensity_var,200)),1);
	rf_i.Set(intensity_mean,sqrt(intensity_var),1);
}

bool MCMC_Intensity_Gaussian::SynthesizeSampling(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = rf_i.Sampling();
	vc_a(1) = vc_a(0);
	vc_a(2) = vc_a(0);

	for (int i=0; i<vc_a.size(); i++)
	{
		if (vc_a(i)<0) vc_a(i) = 0;
		if (vc_a(i)>255) vc_a(i) = 255;
	}

	return true; // in RGB space
}

bool MCMC_Intensity_Gaussian::Reconstruct(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = intensity_mean;
	vc_a(1) = vc_a(0);
	vc_a(2) = vc_a(0);

	return true; // in RGB space
}

mcmcVector<double> MCMC_Intensity_Gaussian::GetAMean(void)
{
	mcmcVector<double> mx_t1(1);

	mx_t1(0) = intensity_mean;

	for (int i=0; i<mx_t1.size(); i++)
	{
		if (mx_t1(i)<0) mx_t1(i) = 0;
		if (mx_t1(i)>255) mx_t1(i) = 255;
	}
	return mx_t1;
}

// get the log of the posterior probability at (row,col)
double MCMC_Intensity_Gaussian::GetLogPostProb(const int row, const int col, mcmcMatrix<int> *pmx_occupied)
{
	assert(prespFilter!=NULL);

	int		i,j;
	double	dis=0.0,count=0.0;

	for (j=row-1; j<=row; j++)
		for (i=col-1; i<=col; i++)
			if (prespFilter->Valid(j,i) && (pmx_occupied==NULL || (pmx_occupied!=NULL && pmx_occupied->Data(j,i)==0)))
			{
				dis += bin_energy.GetBin(prespFilter->ImageRegionLayerIntensity(j,i));
				count = count+1.0;
			}
	if (count<2.0 && pmx_occupied!=NULL)
	{
		dis = 0.0;
		count = 0.0;
		// pixels may be claimed by the curve layer
		for (j=row-1; j<=row+1; j++)
			for (i=col-1; i<=col+1; i++)
				if (prespFilter->Valid(j,i) && pmx_occupied->Data(j,i)==0)
				{
					dis += bin_energy.GetBin(prespFilter->ImageRegionLayerIntensity(j,i));
					count = count+1.0;
				}
	}

	if (count>0)
		dis = dis/count;

	return dis+Complexity();
}


void MCMC_Intensity_Gaussian::OutputInfo(char *info)
{
	if (info!=NULL)
		sprintf(info, "%.2f(%.5f) intensity %.1f:%.1f", log_pos_prob,log_lik/nsize,intensity_mean,intensity_var);
}



//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
// construction & destruction
MCMC_Intensity_Delta::MCMC_Intensity_Delta()
:MCMC_RegionModel()
{
	complexity_prior	= COMPLEXTY_DELTA;
	type				= MODEL_INTENSITY_DELTA;
}

MCMC_Intensity_Delta::~MCMC_Intensity_Delta()
{
	Free();
}

void MCMC_Intensity_Delta::Free(void)
{
	MCMC_RegionModel::Free();
}


// operators
MCMC_Intensity_Delta& MCMC_Intensity_Delta::operator=(const MCMC_Intensity_Delta &x)
{
	if (this!=&x)
	{
		*((MCMC_RegionModel*) this) = *((MCMC_RegionModel*) &x);

		intensity_mean		= x.intensity_mean;
		bin_energy			= x.bin_energy;
	}
	return *this;
}

int MCMC_Intensity_Delta::TotalCodingBits(void)
{
	int bits;

	bits =	bin_energy.TotalCodingBits();	// bins
	return bits;
}


MCMC_Model* MCMC_Intensity_Delta::CopyTo(void)
{
	MCMC_Intensity_Delta *preturn;

	preturn = new MCMC_Intensity_Delta();
	(*preturn) = (*this);

	return ((MCMC_Model*)preturn);
}

void MCMC_Intensity_Delta::SaveToFile(ofstream &strm)
{
	MCMC_RegionModel::SaveToFile(strm);
	strm.write((char*)&intensity_mean, sizeof(intensity_mean));
}

void MCMC_Intensity_Delta::ReadFromFile(ifstream &strm)
{
	MCMC_RegionModel::ReadFromFile(strm);
	strm.read((char*)&intensity_mean, sizeof(intensity_mean));
}

bool MCMC_Intensity_Delta::MeanFromLabel(const int nlabel)
{
	int 	i,j;

	nsize = 0;
	intensity_mean = 0.0;
	bin_energy.Set(INTENSITY_LOW,INTENSITY_HIGH,INTENSITY_DELTA);
	bin_energy.Init(0.0);

	// compute the mean value
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				bin_energy.AddBin(prespFilter->ImageIntensity(j,i),1.0);
				nsize++;

				intensity_mean += prespFilter->ImageIntensity(j,i);
			}
		}
	if (nsize>0)
		intensity_mean = intensity_mean/((double)nsize);
	bin_energy.Normalize();
	for (i=0; i<bin_energy.count; i++)
	{
		if (bin_energy.pbins[i]>1e-20)
			bin_energy.pbins[i] = -log(bin_energy.pbins[i]);
		else
			bin_energy.pbins[i] = 20;
	}

	if (nsize>0)	return true;
	else			return false;
}

bool MCMC_Intensity_Delta::MeanFromPixelSet(void)
{
	int		i,j;

	// initialization
	nsize = 0;
	intensity_mean = 0.0;
	bin_energy.Set(INTENSITY_LOW, INTENSITY_HIGH, INTENSITY_DELTA);
	bin_energy.Init(0.0);

	// get each pixel
	assert(plist_domain != NULL);

	// compute the mean value
	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					bin_energy.AddBin(prespFilter->ImageIntensity(j,i),1.0);
					nsize++;
					intensity_mean += prespFilter->ImageIntensity(j,i);
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (nsize > 0)
		intensity_mean = intensity_mean/((double)nsize);
	bin_energy.Normalize();
	for (i=0; i<bin_energy.count; i++)
	{
		if (bin_energy.pbins[i]>1e-20)
			bin_energy.pbins[i] = -log(bin_energy.pbins[i]);
		else
			bin_energy.pbins[i] = 20;
	}
	
	if (nsize>0)	return true;
	else			return false;
}



// log likelyhood
bool MCMC_Intensity_Delta::LogLikelihoodFromLabel(const int nlabel)
{
	int 	i,j;
	double		dis_uc,din;

	// initialization
	nsize = 0;
	vc_vote.InitValue(0);
	crying_voice = 0.0;
	log_pos_prob = 0.0;

	// the real size of the region
	double	real_size=0.0;
	double	dtemp;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && (*pmx_label)(j,i)==nlabel)
			{
				real_size += 1.0;

				if (ConsiderThisPixel(j,i))
				{
					din = prespFilter->ImageRegionLayerIntensity(j,i);

					dtemp = bin_energy.GetBin(din);

					log_pos_prob += dtemp;
					
					pparticle->AddClusterVote(vc_vote,j,i);

					// compute the crying voice of the regions
					dis_uc = RobustCryingVoice(bin_energy.GetBin(din));
					crying_voice = crying_voice+dis_uc;

					nsize++;
				}
			}
		}
	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	log_lik = log_pos_prob;
	
	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;
	log_pos_prob += ((double)nsize)*Complexity();

	crying_voice = pow(crying_voice,0.6);

	if (nsize>0)	return true;
	else			return false;
}

bool MCMC_Intensity_Delta::LogLikelihoodFromPixelSet(void)
{
	int		i,j;
	double	dis_uc,din;

	// initialization
	nsize = 0;
	vc_vote.InitValue(0);
	crying_voice = 0.0;
	log_pos_prob = 0.0;

	// the real size of the region
	double	real_size,dtemp;

	real_size = 0.0;
	
	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			ppixel_set->energy = 0.0;

			// get each pixel
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				real_size += 1.0;

				if (ConsiderThisPixel(j,i))
				{
					din = prespFilter->ImageRegionLayerIntensity(j,i);

					dtemp = bin_energy.GetBin(din);

					ppixel_set->energy += dtemp;

					log_pos_prob += dtemp;

					// compute the crying voice of the regions
					dis_uc = RobustCryingVoice(bin_energy.GetBin(din));
					crying_voice += dis_uc;

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;
	log_pos_prob += pow((double)nsize, COMPLEXITY_POW)*scale_factor;
	log_pos_prob += ((double)nsize)*Complexity();
	log_lik = log_pos_prob;
	crying_voice = pow(crying_voice, 0.6);

	if (nsize>0)	return true;
	else			return false;
}



void MCMC_Intensity_Delta::UpdateLogLikelihood(const int nlabel, MCMC_PixelSet &list_pixels_old, MCMC_PixelSet &list_pixels_new)
{
	int 	i,j;
	double		din,dtemp,dcomplexity_temp;
	bool		bgo;
	int			k=0;

	if (list_pixels_old.Num()<=0 && list_pixels_new.Num()<=0)
		return;

	log_pos_prob -= pow((double)nsize,COMPLEXITY_POW)*scale_factor;

	dcomplexity_temp = Complexity();
	// delete the old pixels
	bgo = list_pixels_old.GetFirst(j, i);
	while (bgo)
	{
		if (pmx_label->Data(j,i) == nlabel)
		{
			// LogLikelihood in intensity
			din = prespFilter->ImageRegionLayerIntensity(j,i);
			dtemp = bin_energy.GetBin(din);

			log_pos_prob += dtemp+dcomplexity_temp;
			log_lik += dtemp;

			nsize++;
		}
		bgo = list_pixels_old.GetNext(j, i);
	}

	bgo = list_pixels_new.GetFirst(j, i);
	while (bgo)
	{
		if (pmx_label->Data(j,i) == nlabel)
		{
			// LogLikelihood in intensity
			din = prespFilter->ImageRegionLayerIntensity(j,i);
			dtemp = bin_energy.GetBin(din);

			log_pos_prob = log_pos_prob-(dtemp+dcomplexity_temp);
			log_lik -= dtemp;

			nsize--;
		}
		bgo = list_pixels_new.GetNext(j, i);
	}

	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;
}


void MCMC_Intensity_Delta::PrepreForSynthesize(void)
{
	double	value;
	int		i;

	rf_i.Set(INTENSITY_LOW,INTENSITY_HIGH,INTENSITY_DELTA);
	for (i=0; i<bin_energy.count; i++)
	{
		value = rf_i.GetValue(i);
		rf_i.SetOneBin(value, exp(-bin_energy.pbins[i]));
	}
	rf_i.SetCdf();
}

bool MCMC_Intensity_Delta::SynthesizeSampling(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = rf_i.Sampling();
	vc_a(1) = vc_a(0);
	vc_a(2) = vc_a(0);

	return true;		// in RGB space
}

bool MCMC_Intensity_Delta::Reconstruct(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = intensity_mean;
	vc_a(1) = vc_a(0);
	vc_a(2) = vc_a(0);

	return true;		// in RGB space
}

// get the log of the posterior probability at (row,col)
double MCMC_Intensity_Delta::GetLogPostProb(const int row, const int col, mcmcMatrix<int> *pmx_occupied)
{
	assert(prespFilter!=NULL);

	int		i,j;
	double	dis=0.0,count=0.0;

	for (j=row-1; j<=row; j++)
		for (i=col-1; i<=col; i++)
			if (prespFilter->Valid(j,i) &&(pmx_occupied==NULL || (pmx_occupied!=NULL && pmx_occupied->Data(j,i)==0)))
			{
				dis += bin_energy.GetBin(prespFilter->ImageRegionLayerIntensity(j,i));

				count = count+1.0;
			}
	if (count<2.0 && pmx_occupied!=NULL)
	{
		dis = 0.0;
		count = 0.0;

		for (j=row-1; j<=row+1; j++)
			for (i=col-1; i<=col+1; i++)
				if (prespFilter->Valid(j,i) && pmx_occupied->Data(j,i)==0)
				{
					dis += bin_energy.GetBin(prespFilter->ImageRegionLayerIntensity(j,i));

					count = count+1.0;
				}
	}

	if (count>0)
		dis = dis/count;

	return dis+Complexity();
}

void MCMC_Intensity_Delta::OutputInfo(char *info)
{
	if (info!=NULL)
		sprintf(info, "delta %.2f(%.5f)-%.2f", log_pos_prob, log_lik/nsize,intensity_mean);
}

double MCMC_Intensity_Delta::GetPdf(void)
{
	return 1.0;
}

void MCMC_Intensity_Delta::GetDirectValue(mcmcVector<double> &vc, const int row_pixel,const int col_pixel)
{
	vc(0) = prespFilter->ImageIntensity(row_pixel,col_pixel);
}

double MCMC_Intensity_Delta::GetDirectDis(mcmcVector<double> &vc_pixel_intensity)
{
	return pow(vc_pixel_intensity(0)-intensity_mean,2.0);
}




//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
MCMC_Intensity_Surface::MCMC_Intensity_Surface()
:MCMC_RegionModel()
{
	complexity_prior = COMPLEXTY_SURFACE;

	type = MODEL_INTENSITY_SURFACE;

}

MCMC_Intensity_Surface::~MCMC_Intensity_Surface()
{
	Free();
}

void MCMC_Intensity_Surface::Free(void)
{
}

// operators
MCMC_Intensity_Surface& MCMC_Intensity_Surface::operator=(const MCMC_Intensity_Surface &x)
{
	if (this!=&x)
	{
		*((MCMC_RegionModel*) this) = *((MCMC_RegionModel*) &x);

		surface			= x.surface;
		intensity_mean	= x.intensity_mean;
		intensity_var	= x.intensity_var;
		bin_energy		= x.bin_energy;
	}
	return *this;
}

int MCMC_Intensity_Surface::TotalCodingBits(void)
{
	int bits;

	bits = surface.TotalCodingBits();		// control points
	return bits;
}


void MCMC_Intensity_Surface::SaveToFile(ofstream &strm)
{
	MCMC_RegionModel::SaveToFile(strm);

	strm.write((char*)&intensity_mean, sizeof(intensity_mean));
	strm.write((char*)&intensity_var, sizeof(intensity_var));

	double Z=0.0;
	strm.write((char*)&Z, sizeof(Z));
}

void MCMC_Intensity_Surface::ReadFromFile(ifstream &strm)
{
	MCMC_RegionModel::ReadFromFile(strm);
	strm.read((char*)&intensity_mean, sizeof(intensity_mean));
	strm.read((char*)&intensity_var, sizeof(intensity_var));

	double Z;
	strm.read((char*)&Z, sizeof(Z));
}

MCMC_Model* MCMC_Intensity_Surface::CopyTo(void)
{
	MCMC_Intensity_Surface *preturn;

	preturn = new MCMC_Intensity_Surface();
	(*preturn) = (*this);

	return ((MCMC_Model*)preturn);
}



// compute mean and variance for the regions

bool MCMC_Intensity_Surface::MeanFromLabel(const int nlabel)
{

	// compute the intensity means
	int 	i,j;

	// compute the surface first
	surface.Initialize(*pmx_label,nlabel, pmx_taken, prespFilter->ras_intensity);
	surface.Estimate();

	intensity_mean = 0.0;
	nsize = 0;

	// compute the major direction
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				// get intensity
				intensity_mean += prespFilter->ImageIntensity(j,i);

				nsize++;
			}
		}

	if (nsize > 0)
		intensity_mean = intensity_mean/((double)nsize);


	if (nsize>0)	return true;
	else			return false;
}


bool MCMC_Intensity_Surface::MeanFromPixelSet(void)
{
	// compute the surface first
	assert(plist_domain != NULL);

	surface.Initialize(*plist_domain, pmx_taken, prespFilter->ras_intensity);
	surface.Estimate();

	intensity_mean = 0.0;
	nsize = 0;

	int				i,j;
	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// get intensity
					intensity_mean += prespFilter->ImageIntensity(j,i);

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}
	if (nsize > 0)
		intensity_mean = intensity_mean/((double)nsize);

	if (nsize>0)	return true;
	else			return false;
}

bool MCMC_Intensity_Surface::VarianceFromLabel(const int nlabel)
{
	int 	i,j;
	double		din,num_small,small_var,num_big,num_count,intensity_onsurface;

	// initialization
	nsize = 0;
	intensity_var = 0.0;

	num_small = 0.0;
	small_var = 0.0;
	num_big = 0.0;
	num_count = 0.0;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				// projection
				din = prespFilter->ImageIntensity(j,i);
				intensity_onsurface = surface.SurfaceValue(j,i);

				// for robust variance
				if (fabs(din-intensity_onsurface)<5.0)
					num_small += 1.0;
				if (fabs(din-intensity_onsurface)>35.0)
					num_big += 1.0;
				else
				{
					small_var += pow(din-intensity_onsurface,2.0);
					num_count += 1.0;
				}

				// variance of intensity
				intensity_var += pow(din-intensity_onsurface,2.0);

				nsize++;
			}
		}

	if (nsize > 1)
	{
		if (nsize>1000 && num_small/((double)nsize)>0.90 && num_big/((double)nsize)<0.04)
			intensity_var = small_var/num_count;
		else
			intensity_var = intensity_var/((double)nsize-1);
	}
	else
		intensity_var = 1.0;

	if (intensity_var<5.0)
		intensity_var = 5.0;

	// compute quantized energy
	Histogram	hist_temp;
	double		value=0.0,prob,dis;

	bin_energy.Set(0.0,INTENSITY_HIGH,INTENSITY_DELTA);
	hist_temp.Set(0.0,INTENSITY_HIGH,INTENSITY_DELTA);

	while (value<INTENSITY_HIGH+INTENSITY_DELTA)
	{
		dis = pow(value,2.0)/intensity_var/2.0;
		prob = ProbFromEnergy(RobustLoglik(dis)+log(sqrt(2.0*PI))+log(intensity_var)*0.5);
		hist_temp.AddBin(value,prob);
		value += 1.0;
	}
	// hist_temp should not be normlized
	for (i=0; i<hist_temp.count; i++)
	{
		if (hist_temp.pbins[i]>0.0)
			bin_energy.pbins[i] = -log(hist_temp.pbins[i]);
		else
			bin_energy.pbins[i] = 60;
	}

	return true;
}

bool MCMC_Intensity_Surface::VarianceFromPixelSet(void)
{
	double		din,num_small,small_var,num_big,num_count,intensity_onsurface;


	// initialization
	nsize = 0;
	intensity_var = 0.0;
	num_small = 0.0;
	small_var = 0.0;
	num_big = 0.0;
	num_count = 0.0;
	
	int				i,j;
	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// projection
					din = prespFilter->ImageIntensity(j,i);
					intensity_onsurface = surface.SurfaceValue(j,i);

					// for robust variance
					if (fabs(din-intensity_onsurface)<5.0)
						num_small += 1.0;
					if (fabs(din-intensity_onsurface)>35.0)
						num_big += 1.0;
					else
					{
						small_var += pow(din-intensity_onsurface,2.0);
						num_count += 1.0;
					}

					// variance of intensity
					intensity_var += pow(din-intensity_onsurface,2.0);

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (nsize > 1)
	{
		if (nsize>1000 && num_small/((double)nsize)>0.90 && num_big/((double)nsize)<0.04)
			intensity_var = small_var/num_count;
		else
			intensity_var = intensity_var/((double)nsize-1);
	}
	else
		intensity_var = 1.0;

	if (intensity_var < 5.0)
		intensity_var = 5.0;
	else if (intensity_var > 1e5)
		intensity_var = 1e5;

	// compute quantized energy
	Histogram	hist_temp;
	double		value=0.0,prob,dis;

	bin_energy.Set(0.0,INTENSITY_HIGH,INTENSITY_DELTA);
	hist_temp.Set(0.0,INTENSITY_HIGH,INTENSITY_DELTA);

	while (value<INTENSITY_HIGH+INTENSITY_DELTA)
	{
		dis = pow(value,2.0)/intensity_var/2.0;
		prob = ProbFromEnergy(RobustLoglik(dis)+log(sqrt(2.0*PI))+log(intensity_var)*0.5);
		hist_temp.AddBin(value,prob);
		value += 1.0;
	}
	// hist_temp should not be normlized
	for (i=0; i<hist_temp.count; i++)
	{
		if (hist_temp.pbins[i]>0.0)
			bin_energy.pbins[i] = -log(hist_temp.pbins[i]);
		else
			bin_energy.pbins[i] = 60;
	}

	return true;
}

bool MCMC_Intensity_Surface::LogLikelihoodFromLabel(const int nlabel)
{
	int 	i,j;
	double		dis_uc,dis,din,dshift;

	nsize	= 0;
	log_pos_prob = 0.0;

	vc_vote.InitValue(0);
	crying_voice = 0.0;

	double	real_size=0.0;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && (*pmx_label)(j,i)==nlabel)
			{
				real_size += 1.0;

				if (ConsiderThisPixel(j,i))
				{
					// LogLikelihood in intensity
					din = prespFilter->ImageRegionLayerIntensity(j,i);
					dshift = fabs(din-surface.SurfaceValue(j,i));

					dis = pow(dshift,2.0)/intensity_var/2.0;
					log_pos_prob += bin_energy.GetBin(dshift);
					
				
					// compute the crying voice of this region
					dis_uc = RobustCryingVoice(bin_energy.GetBin(dshift)+Complexity());
					crying_voice += dis_uc;

					nsize++;
				}
			}
		}

	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	// the third direction is treated as an independent variable
	log_lik = log_pos_prob;

	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;
	log_pos_prob += ((double)nsize)*Complexity();

	crying_voice = pow(crying_voice,0.6);

	return true;
}


bool MCMC_Intensity_Surface::LogLikelihoodFromPixelSet(void)
{
	double		dis_uc,dis,din,dshift;
	
	// initialization
	nsize	= 0;
	log_pos_prob = 0.0;
	vc_vote.InitValue(0);
	crying_voice = 0.0;

	double	real_size,dtemp;
	int		i,j;
	
	real_size = 0.0;
	
	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			ppixel_set->energy = 0.0;

			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				real_size += 1.0;

				if (ConsiderThisPixel(j,i))
				{
					// LogLikelihood in intensity
					din = prespFilter->ImageRegionLayerIntensity(j,i);
					dshift = fabs(din-surface.SurfaceValue(j,i));

					dis = pow(dshift,2.0)/intensity_var/2.0;
					dtemp = bin_energy.GetBin(dshift);
					log_pos_prob += dtemp;
					ppixel_set->energy += dtemp;

					// compute the crying voice of this region
					dis_uc = RobustCryingVoice(bin_energy.GetBin(dshift)+Complexity());
					crying_voice += dis_uc;

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	// the third direction is treated as an independent variable
	log_lik = log_pos_prob;
	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;
	log_pos_prob += ((double)nsize)*Complexity();
	crying_voice = pow(crying_voice, 0.6);

	return true;
}


void MCMC_Intensity_Surface::UpdateLogLikelihood(const int nlabel, MCMC_PixelSet &list_pixels_old, MCMC_PixelSet &list_pixels_new)
{
	int 	i,j;
	double		dshift,din,dtemp,dcomplexity_temp;
	bool		bgo;

	if (list_pixels_old.Num()<=0 && list_pixels_new.Num()<=0)
		return;

	log_pos_prob -= pow((double)nsize,COMPLEXITY_POW)*scale_factor;

	dcomplexity_temp = Complexity();
	// delete the old pixels
	bgo = list_pixels_old.GetFirst(j, i);
	while (bgo)
	{
		if (pmx_label->Data(j,i) == nlabel)
		{
			// LogLikelihood in intensity
			din = prespFilter->ImageRegionLayerIntensity(j,i);
			dshift = fabs(din-surface.SurfaceValue(j,i));
			dtemp = bin_energy.GetBin(dshift);

			log_pos_prob += dtemp+dcomplexity_temp;
			log_lik += dtemp;

			nsize++;
		}
		bgo = list_pixels_old.GetNext(j, i);
	}

	bgo = list_pixels_new.GetFirst(j, i);
	while (bgo)
	{
		if (pmx_label->Data(j,i) == nlabel)
		{
			// LogLikelihood in intensity
			din = prespFilter->ImageRegionLayerIntensity(j,i);
			dshift = fabs(din-surface.SurfaceValue(j,i));
			dtemp = bin_energy.GetBin(dshift);

			log_pos_prob = log_pos_prob-(dtemp+dcomplexity_temp);
			log_lik -= dtemp;

			nsize--;
		}
		bgo = list_pixels_new.GetNext(j, i);
	}

	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;
}


void MCMC_Intensity_Surface::GetDirectValue(mcmcVector<double> &vc, const int row_pixel,const int col_pixel)
{
	vc(0) = prespFilter->ImageIntensity(row_pixel,col_pixel);
}

double MCMC_Intensity_Surface::GetDirectDis(mcmcVector<double> &vc_intensity)
{
	double	din;

	din = vc_intensity(0);

	return (intensity_mean-din)*(intensity_mean-din); // give lower prevelige over intensity
}


double MCMC_Intensity_Surface::GetPdf(void)
{
	return 1.0;
}

void MCMC_Intensity_Surface::PrepreForSynthesize(void)
{
	rf_i.Set(0.0,sqrt(intensity_var),1);
}

bool MCMC_Intensity_Surface::SynthesizeSampling(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a.SetRowVector(3);

	vc_a(0) = surface.SurfaceValue(row,col)+rf_i.Sampling();
	vc_a(1) = vc_a(0);
	vc_a(2) = vc_a(0);

	for (int i=0; i<vc_a.size(); i++)
	{
		if (vc_a(i)<0) vc_a(i) = 0;
		if (vc_a(i)>255) vc_a(i) = 255;
	}

	return true; // in RGB space
}

bool MCMC_Intensity_Surface::Reconstruct(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = surface.SurfaceValue(row,col);
	vc_a(1) = vc_a(0);
	vc_a(2) = vc_a(0);

	return true;		// in RGB space
}

// get the log of the posterior probability at (row,col)
double MCMC_Intensity_Surface::GetLogPostProb(const int row, const int col, mcmcMatrix<int> *pmx_occupied)
{
	assert(prespFilter!=NULL);

	int		i,j;
	double	dis=0.0,count=0.0,dshift;

	for (j=row-1; j<=row; j++)
		for (i=col-1; i<=col; i++)
			if (prespFilter->Valid(j,i) && (pmx_occupied==NULL || (pmx_occupied!=NULL && pmx_occupied->Data(j,i)==0)))
			{
				dshift = fabs(prespFilter->ImageRegionLayerIntensity(j,i)-surface.SurfaceValue(row,col));

				dis += bin_energy.GetBin(dshift);
				count = count+1.0;
			}
	if (count < 2.0 && pmx_occupied != NULL)
	{
		dis = 0.0;
		count = 0.0;

		for (j=row-1; j<=row+1; j++)
			for (i=col-1; i<=col+1; i++)
				if (prespFilter->Valid(j,i) && pmx_occupied->Data(j,i)==0)
				{
					dshift = fabs(prespFilter->ImageRegionLayerIntensity(j,i)-surface.SurfaceValue(row,col));

					dis += bin_energy.GetBin(dshift);
					count = count+1.0;
				}
	}

	if (count>0)
		dis = dis/count;

	return dis+Complexity();
}


void MCMC_Intensity_Surface::OutputInfo(char *info)
{
	if (info!=NULL)
		sprintf(info, "%.2f(%.5f) surface %.1f:%.1f", log_pos_prob,log_lik/nsize,intensity_mean,intensity_var);
}


////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
//        Models for regions in color images


MCMC_LuvColor_Gaussian::MCMC_LuvColor_Gaussian()
:MCMC_RegionModel()
{
	complexity_prior = COMPLEXTY_SINGLEGAU;

	type = MODEL_LUVCOLOR_GAUSS;

}

MCMC_LuvColor_Gaussian::~MCMC_LuvColor_Gaussian()
{
	Free();
}

void MCMC_LuvColor_Gaussian::Free(void)
{
}

// operators
MCMC_LuvColor_Gaussian& MCMC_LuvColor_Gaussian::operator=(const MCMC_LuvColor_Gaussian &x)
{
	if (this!=&x)
	{
		*((MCMC_RegionModel*) this) = *((MCMC_RegionModel*) &x);

		vc_mean	= x.vc_mean;
		mx_variance	= x.mx_variance;
		
		//compute inverse of the variance matrix
		if ( mx_variance.size()>0 )
		{
			mx_variance_inv = mx_variance.Inv();
			d_std = mx_variance.DetNoZero();
		}
	}
	return *this;
}

int MCMC_LuvColor_Gaussian::TotalCodingBits(void)
{
	int bits;

	bits = sizeof(uchar)* vc_mean.size() +		// mean
		sizeof(float) * mx_variance.rows() * mx_variance.cols();		// variance

	return bits;
}

void MCMC_LuvColor_Gaussian::SaveToFile(ofstream &strm)
{
	MCMC_RegionModel::SaveToFile(strm);

	vc_mean>>strm;
	mx_variance>>strm;
	mx_variance_inv>>strm;

	strm.write((char*)&d_std, sizeof(d_std));
}

void MCMC_LuvColor_Gaussian::ReadFromFile(ifstream &strm)
{
	MCMC_RegionModel::ReadFromFile(strm);

	vc_mean<<strm;
	mx_variance<<strm;
	mx_variance_inv<<strm;

	strm.read((char*)&d_std, sizeof(d_std));
}

MCMC_Model* MCMC_LuvColor_Gaussian::CopyTo(void)
{
	MCMC_LuvColor_Gaussian *preturn;

	preturn = new MCMC_LuvColor_Gaussian();
	(*preturn) = (*this);	

	return ((MCMC_Model*)preturn);
}

// compute mean and variance for the regions
bool MCMC_LuvColor_Gaussian::MeanFromLabel(const int nlabel)
{
	// compute the intensity means
	int 	i,j;

	vc_mean.SetRowVector(3);
	vc_mean.InitValue(0.0);
	nsize = 0;

	// compute the major direction
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				// get L, U, V
				vc_mean(0) += (double)prespFilter->image_original.L(j,i);
				vc_mean(1) += (double)prespFilter->image_original.U(j,i);
				vc_mean(2) += (double)prespFilter->image_original.V(j,i);

				nsize++;
			}
		}

		if (nsize > 0)
		{
			vc_mean = vc_mean / (double)nsize;
		}
			

		// robust mean
		mcmcVector<double> mean_temp;
		mean_temp.SetRowVector(3);
		mean_temp.InitValue(0.0);
		double size_temp1=0.0;
		double size_temp2=0.0;
		double size_temp3=0.0;

		for (j=bound.top; j<=bound.bottom; j++)
			for (i=bound.left; i<=bound.right; i++)
			{
				if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
				{
					// get L,U,V
					if (fabs(prespFilter->image_original.L(j,i)-vc_mean(0))<5.0)
					{
						mean_temp(0) += prespFilter->image_original.L(j,i);						
						size_temp1 += 1.0;
					}

					if (fabs(prespFilter->image_original.U(j,i)-vc_mean(1))<5.0)
					{
						mean_temp(1) += prespFilter->image_original.U(j,i);						
						size_temp2 += 1.0;
					}

					if (fabs(prespFilter->image_original.V(j,i)-vc_mean(2))<5.0)
					{
						mean_temp(2) += prespFilter->image_original.V(j,i);						
						size_temp3 += 1.0;
					}
				}
			}
			if ( nsize>300 )
			{
				if ( size_temp1/nsize>0.90)
				{
					vc_mean(0) = mean_temp(0) / size_temp1;
				}

				if ( size_temp2/nsize>0.90)
				{
					vc_mean(1) = mean_temp(1) / size_temp2;
				}

				if ( size_temp3/nsize>0.90)
				{
					vc_mean(2) = mean_temp(2) / size_temp3;
				}
			}
				

			if (nsize>0)	return true;
			else			return false;
}

bool MCMC_LuvColor_Gaussian::MeanFromPixelSet(void)
{
	// compute the intensity means
	int		i,j;
	bool	bgo,bgo_domain;
	PMCMC_PixelSet ppixel_set;

	vc_mean.SetRowVector(3);
	vc_mean.InitValue(0.0);
	nsize			= 0;

	assert(plist_domain != NULL);

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// get L, U, V
					vc_mean(0) += (double)prespFilter->image_original.L(j,i);
					vc_mean(1) += (double)prespFilter->image_original.U(j,i);
					vc_mean(2) += (double)prespFilter->image_original.V(j,i);
					
					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (nsize > 0)
	{
		vc_mean = vc_mean / (double)nsize;
	}

	// robust mean
	mcmcVector<double> mean_temp;
	mean_temp.SetRowVector(3);
	mean_temp.InitValue(0.0);
	double size_temp1=0.0;
	double size_temp2=0.0;
	double size_temp3=0.0;

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// get L,U,V
					if (fabs(prespFilter->image_original.L(j,i)-vc_mean(0))<5.0)
					{
						mean_temp(0) += prespFilter->image_original.L(j,i);						
						size_temp1 += 1.0;
					}

					if (fabs(prespFilter->image_original.U(j,i)-vc_mean(1))<5.0)
					{
						mean_temp(1) += prespFilter->image_original.U(j,i);						
						size_temp2 += 1.0;
					}

					if (fabs(prespFilter->image_original.V(j,i)-vc_mean(2))<5.0)
					{
						mean_temp(2) += prespFilter->image_original.V(j,i);						
						size_temp3 += 1.0;
					}
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if ( nsize>300 )
	{
		if ( size_temp1/nsize>0.90)
		{
			vc_mean(0) = mean_temp(0) / size_temp1;
		}

		if ( size_temp2/nsize>0.90)
		{
			vc_mean(1) = mean_temp(1) / size_temp2;
		}

		if ( size_temp3/nsize>0.90)
		{
			vc_mean(2) = mean_temp(2) / size_temp3;
		}
	}
	if (nsize>0)	return true;
	else			return false;
}

bool MCMC_LuvColor_Gaussian::VarianceFromLabel(const int nlabel)
{
	int 	i,j,n,m;
	double		num_small,num_count;
	mcmcMatrix<double> small_var;

	// initialization
	nsize = 0;	

	mx_variance.SetDimension(3,3);
	mx_variance.InitValue(0.0);

	small_var.SetDimension(3,3);
	small_var.InitValue(0.0);

	num_small = 0.0;	
	num_count = 0.0;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				for (n=0; n<3; n++)
					for (m=0; m<3; m++)
						mx_variance(n, m) += ((double)(prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[n] - vc_mean(n))) *
						((double)(prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[m] - vc_mean(m)));
				nsize++;
			}
		}

		if (nsize>1)
		{
			mx_variance = mx_variance / (double)(nsize-1);
		}
		
		RobustVariance( mx_variance );
		
		 // make sure the determinate of covariance matrix is not zero
		if (mx_variance.Det()==0)
			mx_variance.Identity();
		mx_variance_inv = mx_variance.Inv();
		d_std = mx_variance.DetNoZero();
		

		return true;
}


template<class Tp>
void printMatrix(mcmcMatrix<Tp> imgSrc, string filename) {

	//Jiadi Yang
	ofstream fout;
	fout.open("F:\\workspace2\\DDMCMC\\Testing\\ddmcmc_" + filename);
	for (int i = 0; i < imgSrc.rows(); i++) {
		for (int j = 0; j < imgSrc.cols(); j++) {
			fout << imgSrc.Data(i, j) << " ";
		}
		fout << endl;
	}
	fout.close();
	// 
}

bool MCMC_LuvColor_Gaussian::VarianceFromPixelSet(void)
{
	int		i,j,n,m;
	double	num_small,num_count;
	bool	bgo;

	// initialization
	mcmcMatrix<double> small_var;

	// initialization
	nsize = 0;	

	mx_variance.SetDimension(3,3);
	mx_variance.InitValue(0.0);

	small_var.SetDimension(3,3);
	small_var.InitValue(0.0);

	num_small = 0.0;	
	num_count = 0.0;

	assert(plist_domain != NULL);

	bool			bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					for (n=0; n<3; n++)
					for (m = 0; m < 3; m++) {
						mx_variance(n, m) += ((double)(prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[n] - vc_mean(n))) *
							((double)(prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[m] - vc_mean(m)));
					}
					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (nsize>1)
	{
		mx_variance = mx_variance / (double)(nsize-1);
	}

	RobustVariance( mx_variance );

	// make sure the determinate of covariance matrix is not zero
	if (mx_variance.Det()==0)
		mx_variance.Identity();
	mx_variance_inv = mx_variance.Inv();
	d_std = mx_variance.DetNoZero();

	
	return true;
}


bool MCMC_LuvColor_Gaussian::LogLikelihoodFromLabel(const int nlabel)
{
	int 		i,j,m;
	mcmcMatrix<double>	I_u(3,1);
	double			dis_uc,dis,dtmp[3];
	
	nsize	= 0;
	log_pos_prob = 0.0;

	mx_variance_inv = mx_variance.Inv();
	vc_vote.InitValue(0);
	crying_voice = 0.0;

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				I_u(0) = (double)prespFilter->image_original.L(j,i)-vc_mean(0);
				I_u(1) = (double)prespFilter->image_original.U(j,i)-vc_mean(1);
				I_u(2) = (double)prespFilter->image_original.V(j,i)-vc_mean(2);

				//dis = I_u.T()*mx_variance_inv*I_u/2.0; to speed up
				for (m=0; m<3; m++)
					dtmp[m] = I_u(0)*mx_variance_inv(0,m) + I_u(1)*mx_variance_inv(1,m) +
					I_u(2)*mx_variance_inv(2,m);

				dis = (dtmp[0]*I_u(0)+dtmp[1]*I_u(1)+dtmp[2]*I_u(2))/2.0;

				log_pos_prob += RobustLoglik(dis);
				pparticle->AddClusterVote(vc_vote,j,i);

				// compute the crying voice of this region
				dis_uc = RobustCryingVoice(RobustLoglik(dis));
				crying_voice += dis_uc;

				nsize++;
			}
		}

	
		log_pos_prob += log(sqrt(2.0*PI))*3.0*nsize;

		log_pos_prob += log(mx_variance.DetNoZero())*0.5*nsize;

		// penalize regions with few valid pixels
		if ( nsize <= MINIMAL_SIZE)
			log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

		// the third direction is treated as an independent variable
		log_lik = log_pos_prob;

		log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;

		log_pos_prob += ((double)nsize)*Complexity();

		crying_voice = pow(crying_voice,0.6);

		
		return true;
}

bool MCMC_LuvColor_Gaussian::LogLikelihoodFromPixelSet(void)
{
	int 		i,j,m;
	mcmcMatrix<double>	I_u(3,1);
	double			dis_uc,dis,dtmp[3];

	// initialization
	nsize	= 0;
	log_pos_prob = 0.0;
	vc_vote.InitValue(0);
	crying_voice = 0.0;

	// the real size of the region
	double	real_size;
	real_size = 0.0;

	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			ppixel_set->energy = 0.0;

			// get each pixel
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				real_size += 1.0;

				if (ConsiderThisPixel(j,i))
				{
					
						I_u(0) = (double)prespFilter->image_original.L(j,i)-vc_mean(0);
						I_u(1) = (double)prespFilter->image_original.U(j,i)-vc_mean(1);
						I_u(2) = (double)prespFilter->image_original.V(j,i)-vc_mean(2);

						//dis = I_u.T()*mx_variance_inv*I_u/2.0; to speed up
						for (m=0; m<3; m++)
							dtmp[m] = I_u(0)*mx_variance_inv(0,m) + I_u(1)*mx_variance_inv(1,m) +
							I_u(2)*mx_variance_inv(2,m);

						dis = (dtmp[0]*I_u(0)+dtmp[1]*I_u(1)+dtmp[2]*I_u(2))/2.0;

						log_pos_prob += RobustLoglik(dis);
						pparticle->AddClusterVote(vc_vote,j,i);

						// compute the crying voice of this region
						dis_uc = RobustCryingVoice(RobustLoglik(dis));
						crying_voice += dis_uc;

						nsize++;
					
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	log_pos_prob += log(sqrt(2.0*PI))*3.0*nsize;

	log_pos_prob += log(mx_variance.DetNoZero())*0.5*nsize;


	// penalize regions with few valid pixels
	if (nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	// the third direction is treated as an independent variable
	log_lik = log_pos_prob;
	log_pos_prob += pow((double)nsize, COMPLEXITY_POW)*scale_factor;
	log_pos_prob += ((double)nsize)*Complexity();
	crying_voice = pow(crying_voice, 0.6);

	return true;
}


void MCMC_LuvColor_Gaussian::GetDirectValue(mcmcVector<double> &vc, const int row_pixel,const int col_pixel)
{
	vc(0) = prespFilter->image_original.L(row_pixel,col_pixel);
	vc(1) = prespFilter->image_original.U(row_pixel,col_pixel);
	vc(2) = prespFilter->image_original.V(row_pixel,col_pixel);
}

double MCMC_LuvColor_Gaussian::GetDirectDis(mcmcVector<double> &vc_pixel_luv)
{
	double dtemp=0.0;

	for (int i=0; i<vc_pixel_luv.size(); i++)
		dtemp += pow(vc_pixel_luv(i)-vc_mean(i),2.0);
	return dtemp;
}

void MCMC_LuvColor_Gaussian::DiffuseComputedMean(const int nlabel)
{
	int ipar;

	// compute all the votes
	ComputeVote(nlabel);
	// choose a particle
	ipar = pparticle->ChooseParticle(vc_vote);
	vc_mean	= DiffuseTo(*pparticle,vc_mean,ipar);
}

double MCMC_LuvColor_Gaussian::GetPdf(void)
{	
	assert(pparticle!=NULL);
	
	return pparticle->GetPdf(vc_vote,vc_mean);
}

void MCMC_LuvColor_Gaussian::PrepreForSynthesize(void)
{
	mcmcMatrix<double> mx_eigenvalues;

	mx_variance.Eigen(mx_eigenvalues, mx_eigehvectors);
	rf_l.Set(0.0,sqrt(fabs(mx_eigenvalues(0))),1);
	rf_u.Set(0.0,sqrt(fabs(mx_eigenvalues(1))),1);
	rf_v.Set(0.0,sqrt(fabs(mx_eigenvalues(2))),1);
}

bool MCMC_LuvColor_Gaussian::SynthesizeSampling(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = (int)rf_l.Sampling();
	vc_a(1) = (int)rf_u.Sampling();
	vc_a(2) = (int)rf_v.Sampling();
	mcmcMatrix<double> vc_temp;
	vc_temp.SetDimension(vc_a.cols(), vc_a.rows());
	vc_temp = mx_eigehvectors*(vc_a.T()); 
	CopyMatrix(vc_a,vc_temp.T());


	vc_a(0) += vc_mean(0);
	vc_a(1) += vc_mean(1);
	vc_a(2) += vc_mean(2);

	return false; // in Luv space
}

bool MCMC_LuvColor_Gaussian::Reconstruct(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = vc_mean(0);
	vc_a(1) = vc_mean(1);
	vc_a(2) = vc_mean(2);

	return false;	// in Luv space
}

mcmcVector<double> MCMC_LuvColor_Gaussian::GetAMean(void)
{
	return vc_mean;
}

// get the log of the posterior probability at (row,col)
double MCMC_LuvColor_Gaussian::GetLogPostProb(const int row, const int col, mcmcMatrix<int> *pmx_occupied)
{
	mcmcVector<double>	I_u,dis,vc_pixel(3);
	double			std_temp;

	assert(prespFilter!=NULL);
	
	if ( !prespFilter->Valid(row, col) ) return 0.0; // add by wtf

	vc_pixel(0) = prespFilter->image_original.L(row,col);
	vc_pixel(1) = prespFilter->image_original.U(row,col);
	vc_pixel(2) = prespFilter->image_original.V(row,col); // change by wtf: row to row -1

	I_u = vc_pixel-vc_mean;

	if (d_std >= 100.0)
	{
		dis = I_u*mx_variance_inv*(I_u.T())/2.0;
		std_temp = d_std;
	}
	else
	{
		// if variance is too small then the force will be falsified
		mcmcMatrix<double> mx_temp;
		mx_temp = mx_variance_inv/pow(100.0/d_std,1.0/((double) mx_variance_inv.rows()));
		dis = I_u*mx_temp*(I_u.T())/2.0;
		std_temp = 100.0;
	}

	return	log(sqrt(2.0*PI))*3.0+log(std_temp)*0.5*2.0+dis.Data(0);
}


double MCMC_LuvColor_Gaussian::RobustLoglik(const double dvalue)
{
	if (dvalue>15.0)
		return 15.0;
	else
		return dvalue;
}

void MCMC_LuvColor_Gaussian::OutputInfo(char *info)
{
	if (info!=NULL)
		sprintf(info, "%.2f(%.5f) gaussian %.1f:%.1f:%.1f--%.1f", log_pos_prob,log_lik/nsize,vc_mean(0),
		vc_mean(1),vc_mean(2),d_std);
}



//
MCMC_LuvColor_MixtureGaussians::MCMC_LuvColor_MixtureGaussians()
:MCMC_RegionModel()
{
	complexity_prior = COMPLEXTY_MIXTUREGAU;
	type = MODEL_LUVCOLOR_MIXTUREGAUSS;

}

MCMC_LuvColor_MixtureGaussians::~MCMC_LuvColor_MixtureGaussians()
{
	Free();
}

void MCMC_LuvColor_MixtureGaussians::Free(void)
{
}

// operators
MCMC_LuvColor_MixtureGaussians& MCMC_LuvColor_MixtureGaussians::operator=(const MCMC_LuvColor_MixtureGaussians &x)
{
	if (this!=&x)
	{
		*((MCMC_RegionModel*) this) = *((MCMC_RegionModel*) &x);

		for (int i=0; i<2; i++)
		{
			vc_means[i]		= x.vc_means[i];
			mx_variances[i] = x.mx_variances[i];
			alphs[i]		= x.alphs[i];
			// compute inverse of the variance matrix
			if (mx_variances[i].size()>0)
			{
				mx_variance_invs[i] = mx_variances[i].Inv();
				stds[i] = mx_variances[i].DetNoZero();
			}
		}
	}
	return *this;
}

int MCMC_LuvColor_MixtureGaussians::TotalCodingBits(void)
{
	int bits;

	bits = 0;
	for ( int i=0; i<2; i++ )
	{
		bits += sizeof(uchar)* vc_means[i].size() +		// mean
			sizeof(float) * mx_variances[i].rows() * mx_variances[i].cols();		// variance
	}
	

	return bits;
}

void MCMC_LuvColor_MixtureGaussians::SaveToFile(ofstream &strm)
{
	MCMC_RegionModel::SaveToFile(strm);

	for (int i=0; i<2; i++)
	{
		vc_means[i]>>strm;
		mx_variances[i]>>strm;
		strm.write((char *)&(alphs[i]),sizeof(alphs[i]));
		mx_variance_invs[i]>>strm;
		strm.write((char *)&(stds[i]),sizeof(stds[i]));
	}
}

void MCMC_LuvColor_MixtureGaussians::ReadFromFile(ifstream &strm)
{
	MCMC_RegionModel::ReadFromFile(strm);

	for (int i=0; i<2; i++)
	{
		vc_means[i]<<strm;
		mx_variances[i]<<strm;
		strm.read((char *)&(alphs[i]),sizeof(alphs[i]));
		mx_variance_invs[i]<<strm;
		strm.read((char *)&(stds[i]),sizeof(stds[i]));
	}
}

MCMC_Model* MCMC_LuvColor_MixtureGaussians::CopyTo(void)
{
	MCMC_LuvColor_MixtureGaussians *preturn;

	preturn = new MCMC_LuvColor_MixtureGaussians();
	(*preturn) = (*this);

	return ((MCMC_Model*)preturn);
}

// compute mean and variance for the regions
bool MCMC_LuvColor_MixtureGaussians::MeanFromLabel(const int nlabel)
{	
	int 		i,j;
	int				m;
	double	vc_num_pixel1,vc_num_pixel2;
	double	vc_dtemp1,vc_dtemp2;

	// initialization
	
	for (i=0; i<2; i++)
	{
		vc_means[i].SetRowVector(3);
		vc_means[i].InitValue(0.0);
	}
	
	
	// get two initial values
	double dist;
	vc_dtemp1 = 99999.0;
	vc_dtemp2 = 0.0;

	// get two initial values
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel )
			{			

				dist = pow(prespFilter->image_original.L(j,i)+100.0,2.0)+pow(prespFilter->image_original.U(j,i)+100.0,2.0)+pow(prespFilter->image_original.V(j,i)+100.0,2.0);
				if (dist<vc_dtemp1)
				{
					// the first initial value
					vc_dtemp1 = dist;
					vc_means[0](0) = prespFilter->image_original.L(j,i);
					vc_means[0](1) = prespFilter->image_original.U(j,i);
					vc_means[0](2) = prespFilter->image_original.V(j,i);
				}
				// the second initial value
				if (dist>vc_dtemp2)
				{
					// the second initial value
					vc_dtemp2 = dist;
					vc_means[1](0) = prespFilter->image_original.L(j,i);
					vc_means[1](1) = prespFilter->image_original.U(j,i);
					vc_means[1](2) = prespFilter->image_original.V(j,i);
				}
			}
		}

		// K_Mean clustering
		for (int step=0; step<7; step++)
		{
			vc_num_pixel1 = 0.0;
			vc_num_pixel2 = 0.0;

			for (j=bound.top; j<=bound.bottom; j++)
				for (i=bound.left; i<=bound.right; i++)
				{
					if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel )
					{						

						if (CloserMode(j,i,vc_means[0],vc_means[1])==0)
						{
							for (m=0; m<3; m++)
								vc_means[0](m) = (vc_means[0](m)*(vc_num_pixel1+1) +
													prespFilter->image_original.LUV.at<cv::Vec3i>(j,i)[m])/(vc_num_pixel1+2);
							vc_num_pixel1 += 1.0;
						}
						else
						{
							for (m=0; m<3; m++)
								vc_means[1](m) = (vc_means[1](m)*(vc_num_pixel2+1)+
								prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[m]) / (vc_num_pixel2 + 2);
							vc_num_pixel2 += 1.0;
						}
					}
				}
		}

		return true;
}

bool MCMC_LuvColor_MixtureGaussians::MeanFromPixelSet(void)
{
	// compute the intensity means
	int 		i,j;	
	int  m;
	double	vc_num_pixel1,vc_num_pixel2;
	double	vc_dtemp1,vc_dtemp2;

	bool	bgo,bgo_domain;
	PMCMC_PixelSet ppixel_set;

	for (i=0; i<2; i++)
	{
		vc_means[i].SetRowVector(3);
		vc_means[i].InitValue(0.0);
	}
	

	// get two initial values
	double dist;
	vc_dtemp1 = 99999.0;
	vc_dtemp2 = 0.0;

	assert(plist_domain != NULL);

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					dist = pow(prespFilter->image_original.L(j,i)+100.0,2.0)+pow(prespFilter->image_original.U(j,i)+100.0,2.0)+pow(prespFilter->image_original.V(j,i)+100.0,2.0);
					if (dist<vc_dtemp1)
					{
						// the first initial value
						vc_dtemp1 = dist;
						vc_means[0](0) = prespFilter->image_original.L(j,i);
						vc_means[0](1) = prespFilter->image_original.U(j,i);
						vc_means[0](2) = prespFilter->image_original.V(j,i);
					}
					// the second initial value
					if (dist>vc_dtemp2)
					{
						// the second initial value
						vc_dtemp2 = dist;
						vc_means[1](0) = prespFilter->image_original.L(j,i);
						vc_means[1](1) = prespFilter->image_original.U(j,i);
						vc_means[1](2) = prespFilter->image_original.V(j,i);
					}
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	// K_Mean clustering
	for (int step=0; step<7; step++)
	{
		vc_num_pixel1 = 0.0;
		vc_num_pixel2 = 0.0;

		bgo_domain = plist_domain->GetFirst(ppixel_set);
		while (bgo_domain)
		{
			if (ppixel_set != NULL)
			{
				bgo = ppixel_set->GetFirst(j, i);
				while (bgo)
				{
					if (ConsiderThisPixel(j,i))
					{									

						if (!CloserMode(j,i,vc_means[0],vc_means[1]))
						{
							for (m=0; m<3; m++)
								vc_means[0](m) = (vc_means[0](m)*(vc_num_pixel1+1) +
								prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[m]) / (vc_num_pixel1 + 2);
							vc_num_pixel1 += 1.0;
						}
						else
						{
							for (m=0; m<3; m++)
								vc_means[1](m) = (vc_means[1](m)*(vc_num_pixel2+1)+
								prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[m]) / (vc_num_pixel2 + 2);
							vc_num_pixel2 += 1.0;
						}
						
					}
					bgo = ppixel_set->GetNext(j, i);
				}
			}
			bgo_domain = plist_domain->GetNext(ppixel_set);
		}

	}

	return true;
	
}

bool MCMC_LuvColor_MixtureGaussians::VarianceFromLabel(const int nlabel)
{
	int  i,j,l,n,m;
	int	mx_counts[2];

	// initializations
	nsize = 0;
	mx_counts[0] = 0;
	mx_counts[1] = 0;

	for (i=0; i<2; i++)
	{
		mx_variances[i].SetDimension(3,3);
		mx_variances[i].InitValue(0.0);
	}

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				l = CloserMode(j,i,vc_means[0],vc_means[1]);

				for (n=0; n<3; n++)
					for (m=0; m<3; m++)
						mx_variances[l](n, m) += ((double)(prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[n] - vc_means[l](n))) *
						((double)(prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[m] - vc_means[l](m)));

				mx_counts[l]++;				
			}
		}

	for ( l=0; l<2; l++ )
	{
		if ( mx_counts[l]>1 )
			mx_variances[l] = mx_variances[l] / (mx_counts[l]-1);
		
		RobustVariance(mx_variances[l]);

		alphs[l] = ((double)mx_counts[l])/(mx_counts[0]+mx_counts[1]);
	}


	return true;
}

bool MCMC_LuvColor_MixtureGaussians::VarianceFromPixelSet(void)
{
	int  i,j,l,n,m;
	int	mx_counts[2];

	// initializations
	nsize = 0;
	mx_counts[0] = 0;
	mx_counts[1] = 0;

	for (i=0; i<2; i++)
	{
		mx_variances[i].SetDimension(3,3);
		mx_variances[i].InitValue(0.0);
	}

	assert(plist_domain != NULL);

	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					l = CloserMode(j,i,vc_means[0],vc_means[1]);

					for (n=0; n<3; n++)
						for (m=0; m<3; m++)
							mx_variances[l](n, m) += ((double)(prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[n] - vc_means[l](n))) *
							((double)(prespFilter->image_original.LUV.at<cv::Vec3i>(j, i)[m] - vc_means[l](m)));

					mx_counts[l]++;				
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}
	
	for ( l=0; l<2; l++ )
	{
		if ( mx_counts[l]>1 )
			mx_variances[l] = mx_variances[l] / (mx_counts[l] - 1);

		RobustVariance( mx_variances[l] );

		alphs[l] = ((double)mx_counts[l])/(mx_counts[0]+mx_counts[1]);
	}
    
	return true;
}


bool MCMC_LuvColor_MixtureGaussians::LogLikelihoodFromLabel(const int nlabel)
{
	int 		i,j;
	int				l;	
	mcmcMatrix<double>	*p_mx_variance_invs[2];
	mcmcMatrix<double>	I_u(3,1),dis;
	double			dis_uc,Log_temp[2],eng_temp;
	double          mx_stds[2];


	//model complexity
	log_pos_prob = 0.0;
	crying_voice = 0.0;
	vc_vote.InitValue(0);
	mx_stds[0] = 0;
	mx_stds[1] = 0;

	nsize = 0;

	for ( l=0; l<2; l++ )
	{
		p_mx_variance_invs[l] = new mcmcMatrix<double>;
		*(p_mx_variance_invs[l]) = mx_variances[l].Inv();
		mx_stds[l] = mx_variances[l].DetNoZero();  

	}
	

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{
				for ( l=0; l<2; l++ )
				{
					I_u(0) = (double)prespFilter->image_original.L(j,i);
					I_u(1) = (double)prespFilter->image_original.U(j,i);
					I_u(2) = (double)prespFilter->image_original.V(j,i);

					I_u = I_u - vc_means[l].T();
					dis = I_u.T() * (*p_mx_variance_invs[l]) * I_u / 2.0;

					Log_temp[l] = -log(alphs[l]) + 0.5 * log(mx_stds[l]) + log(sqrt(2.0*PI)) * 3.0 + dis(0);

				}
				
				if ( Log_temp[0]>600 || Log_temp[1]>600 )
				{
					eng_temp = min(Log_temp[0],Log_temp[1]);
				}
				else
				{
					eng_temp = -log(exp(-Log_temp[0])+exp(-Log_temp[1]));
				}

				log_pos_prob += RobustLoglik(eng_temp);

				pparticle->AddClusterVote(vc_vote, j, i);

				// compute the crying voice of this region
				dis_uc = RobustCryingVoice(eng_temp);
				crying_voice += dis_uc;

				nsize++;
			}
		}


		
		// penalize regions with few valid pixels
		if ( nsize <= MINIMAL_SIZE)
			log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

		// the third direction is treated as an independent variable
		log_lik = log_pos_prob;

		log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;

		log_pos_prob += ((double)nsize)*Complexity();

		crying_voice = pow(crying_voice,0.6);

		delete [] p_mx_variance_invs[0];
		delete [] p_mx_variance_invs[1];

		return true;
}

bool MCMC_LuvColor_MixtureGaussians::LogLikelihoodFromPixelSet(void)
{
	int 		i,j;
	int				l;	
	mcmcMatrix<double>	*p_mx_variance_invs[2];
	mcmcMatrix<double>	I_u(3,1);
	double			dis, dis_uc,Log_temp[2],eng_temp;
	double          mx_stds[2];


	//model complexity
	log_pos_prob = 0.0;
	crying_voice = 0.0;
	vc_vote.InitValue(0);
	mx_stds[0] = 0;
	mx_stds[1] = 0;

	nsize = 0;

	for ( l=0; l<2; l++ )
	{
		p_mx_variance_invs[l] = new mcmcMatrix<double>;
		*(p_mx_variance_invs[l]) = mx_variances[l].Inv();
		mx_stds[l] = mx_variances[l].DetNoZero();  

	}

	// the real size of the region
	double	real_size;
	real_size = 0.0;

	bool			bgo,bgo_domain;
	PMCMC_PixelSet	ppixel_set;

	// get each pixel
	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			ppixel_set->energy = 0.0;

			// get each pixel
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				real_size += 1.0;

				if (ConsiderThisPixel(j,i))
				{

					for ( l=0; l<2; l++ )
					{
						I_u(0) = (double)prespFilter->image_original.L(j, i) - vc_means[l](0);
						I_u(1) = (double)prespFilter->image_original.U(j, i) - vc_means[l](1);
						I_u(2) = (double)prespFilter->image_original.V(j, i) - vc_means[l](2);
						dis = (I_u(0)*I_u(0)*(*p_mx_variance_invs[l])(0, 0)
							+ I_u(0)*I_u(1)*((*p_mx_variance_invs[l])(0, 1) + (*p_mx_variance_invs[l])(1, 0))
							+ I_u(1)*I_u(1)*(*p_mx_variance_invs[l])(1, 1)
							+ I_u(1)*I_u(2)*((*p_mx_variance_invs[l])(1, 2) + (*p_mx_variance_invs[l])(2, 1))
							+ I_u(2)*I_u(2)*(*p_mx_variance_invs[l])(2, 2) 
							+ I_u(0)*I_u(2)*((*p_mx_variance_invs[l])(0, 2) + (*p_mx_variance_invs[l])(2, 0))) / 2;


						Log_temp[l] = -log(alphs[l]) + 0.5 * log(mx_stds[l]) + log(sqrt(2.0*PI)) * 3.0 + dis;

					}

					if ( Log_temp[0]>600 || Log_temp[1]>600 )
					{
						eng_temp = min(Log_temp[0],Log_temp[1]);
					}
					else
					{
						eng_temp = -log(exp(-Log_temp[0])+exp(-Log_temp[1]));
					}

					log_pos_prob += RobustLoglik(eng_temp);

					pparticle->AddClusterVote(vc_vote, j, i);

					// compute the crying voice of this region
					dis_uc = RobustCryingVoice(eng_temp);
					crying_voice += dis_uc;

					nsize++;

				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	// penalize regions with few valid pixels
	if ( nsize <= MINIMAL_SIZE)
		log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

	// the third direction is treated as an independent variable
	log_lik = log_pos_prob;

	log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor;

	log_pos_prob += ((double)nsize)*Complexity();

	crying_voice = pow(crying_voice,0.6);

	delete p_mx_variance_invs[0];
	delete p_mx_variance_invs[1];

	return true;
}


void MCMC_LuvColor_MixtureGaussians::GetDirectValue(mcmcVector<double> &vc, const int row_pixel,const int col_pixel)
{
	vc(0) = prespFilter->image_original.L(row_pixel,col_pixel);
	vc(1) = prespFilter->image_original.U(row_pixel,col_pixel);
	vc(2) = prespFilter->image_original.V(row_pixel,col_pixel);
}

double MCMC_LuvColor_MixtureGaussians::GetDirectDis(mcmcVector<double> &vc_pixel_luv)
{
	int k=0;
	double dtemp=0.0;

	k = CloserMode(vc_pixel_luv);

	for (int i=0; i<vc_pixel_luv.size(); i++)
		dtemp += pow(vc_pixel_luv(i)-vc_means[k](i),2.0);
	return dtemp;
}

void MCMC_LuvColor_MixtureGaussians::DiffuseComputedMean(const int nlabel)
{
	int ipar1,ipar2,n=0;

	// compute all the votes
	ComputeVote(nlabel);
	// choose a particle
	ipar1 = pparticle->ChooseParticle(vc_vote);
	ipar2 = ipar1;

	while (ipar2==ipar1 && n<100)
	{
		ipar2 = pparticle->ChooseParticle(vc_vote);
		n++;
	}

	if (pparticle->CloserParticle(vc_means[0],vc_means[1],ipar1)==1)
	{
		// switch
		int temp=ipar2;
		ipar2 = ipar1;
		ipar1 = temp;
	}

	// diffuse to the position inside the particle
	vc_means[0]	= DiffuseTo(*pparticle,vc_means[0],ipar1);	
	vc_means[1]	= DiffuseTo(*pparticle,vc_means[1],ipar2);
	
}

double MCMC_LuvColor_MixtureGaussians::GetPdf(void)
{	
	assert(pparticle!=NULL);

	return pparticle->GetPdf(vc_vote,vc_means[0]) * pparticle->GetPdf(vc_vote,vc_means[1]);
}

void MCMC_LuvColor_MixtureGaussians::PrepreForSynthesize(void)
{
	mcmcMatrix<double> mx_eigenvalues;

	// random field for the two modes
	rf_mode.Set(0,1,1);

	for (int i=0; i<2; i++)
	{
		rf_mode.SetOneBin(i,alphs[i]);

		mx_variances[i].Eigen(mx_eigenvalues, mx_eigehvectors[i]);
		rf_l[i].Set(0.0,sqrt(fabs(mx_eigenvalues(0))),1);
		rf_u[i].Set(0.0,sqrt(fabs(mx_eigenvalues(1))),1);
		rf_v[i].Set(0.0,sqrt(fabs(mx_eigenvalues(2))),1);
	}
	rf_mode.SetCdf();
}

bool MCMC_LuvColor_MixtureGaussians::SynthesizeSampling(mcmcVector<double> &vc_a, const int row, const int col)
{
	int i;

	i = (int)rf_mode.Sampling();

	vc_a(0) = (int)rf_l[i].Sampling();
	vc_a(1) = (int)rf_u[i].Sampling();
	vc_a(2) = (int)rf_v[i].Sampling();

	mcmcMatrix<double> vc_temp;
	vc_temp.SetDimension(vc_a.cols(), vc_a.rows());
	vc_temp = mx_eigehvectors[i]*(vc_a.T()); 
	CopyMatrix(vc_a,vc_temp.T());

	//vc_a = mx_eigehvectors[i]*vc_a;

	vc_a(0) += vc_means[i](0);
	vc_a(1) += vc_means[i](1);
	vc_a(2) += vc_means[i](2);

	return false; // in luv space
}

bool MCMC_LuvColor_MixtureGaussians::Reconstruct(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = alphs[0]*vc_means[0](0)*alphs[1]*vc_means[1](0);
	vc_a(1) = alphs[0]*vc_means[0](1)*alphs[1]*vc_means[1](1);
	vc_a(2) = alphs[0]*vc_means[0](2)*alphs[1]*vc_means[1](2);

	return false;	// in Luv space
}

mcmcVector<double> MCMC_LuvColor_MixtureGaussians::GetAMean(void)
{
	return vc_means[0]*alphs[0]+vc_means[1]*alphs[1];;
}

// get the log of the posterior probability at (row,col)
double MCMC_LuvColor_MixtureGaussians::GetLogPostProb(const int row_pixel, const int col_pixel, mcmcMatrix<int> *pmx_occupied)
{
	mcmcVector<double>	I_u, vc_pixel(3);
	int				i;
	double			dis, dtemp,std_temp,Log_temp[2];

	assert(prespFilter!=NULL);
	
	if ( !prespFilter->Valid(row_pixel, col_pixel) ) return 0.0; // add by wtf

	vc_pixel(0) = prespFilter->image_original.L(row_pixel,col_pixel);
	vc_pixel(1) = prespFilter->image_original.U(row_pixel,col_pixel);
	vc_pixel(2) = prespFilter->image_original.V(row_pixel,col_pixel); //change by wtf: add -1

	for (i=0; i<=1; i++)
	{
		I_u = vc_pixel-vc_means[i];

		if (stds[i] >= 100.0)
		{

			dis = (I_u(0)*I_u(0)*mx_variance_invs[i](0, 0)
				+ I_u(0)*I_u(1)*(mx_variance_invs[i](0, 1) + mx_variance_invs[i](1, 0))
				+ I_u(1)*I_u(1)*mx_variance_invs[i](1, 1)
				+ I_u(1)*I_u(2)*(mx_variance_invs[i](1, 2) + mx_variance_invs[i](2, 1))
				+ I_u(2)*I_u(2)*mx_variance_invs[i](2, 2)
				+ I_u(0)*I_u(2)*(mx_variance_invs[i](0, 2) + mx_variance_invs[i](2, 0))) / 2;
			std_temp = stds[i];
		}
		else
		{
			// if variance is too small then the force will be falsified
			mcmcMatrix<double> mx_temp;
			mx_temp = mx_variance_invs[i]/pow(100.0/stds[i],1.0/((double) mx_variance_invs[i].rows()));
			dis = (I_u(0)*I_u(0)*mx_temp(0, 0)
				+ I_u(0)*I_u(1)*(mx_temp(0, 1) + mx_temp(1, 0))
				+ I_u(1)*I_u(1)*mx_temp(1, 1)
				+ I_u(1)*I_u(2)*(mx_temp(1, 2) + mx_temp(2, 1))
				+ I_u(2)*I_u(2)*mx_temp(2, 2)
				+ I_u(0)*I_u(2)*(mx_temp(0, 2) + mx_temp(2, 0))) / 2;
			std_temp = 100.0;
		}
		Log_temp[i] = log(std_temp)*0.5*2.0+dis-log(alphs[i])+log(sqrt(2.0*PI))*3.0; 
	}

	if (Log_temp[0]>600 || Log_temp[1]>600)
		dtemp = min(Log_temp[0],Log_temp[1]);
	else
		dtemp = -log(exp(-Log_temp[0])+exp(-Log_temp[1]));
	return dtemp;
}


double MCMC_LuvColor_MixtureGaussians::RobustLoglik(const double dvalue)
{
	if (dvalue>15.0)
		return 15.0;
	else
		return dvalue;
}

void MCMC_LuvColor_MixtureGaussians::OutputInfo(char *info)
{
	if (info!=NULL)
		sprintf(info, "%.2f(%.5f) mixture %.1f--%.1f", log_pos_prob,log_lik/nsize,stds[0],stds[1]);
}

int	MCMC_LuvColor_MixtureGaussians::CloserMode(const int j, const int i, mcmcVector<double> &vc_m1, mcmcVector<double> &vc_m2)
{
	double dist1,dist2;

	dist1 = pow(prespFilter->image_original.L(j,i)-vc_m1(0),2.0)+pow(prespFilter->image_original.U(j,i)-vc_m1(1),2.0)+pow(prespFilter->image_original.V(j,i)-vc_m1(2),2.0);
	dist2 = pow(prespFilter->image_original.L(j,i)-vc_m2(0),2.0)+pow(prespFilter->image_original.U(j,i)-vc_m2(1),2.0)+pow(prespFilter->image_original.V(j,i)-vc_m2(2),2.0);
	if (dist1<dist2)
		return 0;
	else
		return 1;
}

int	MCMC_LuvColor_MixtureGaussians::CloserMode(mcmcVector<double> &vc_pixel)
{
	double dist1,dist2;

	dist1 = pow(vc_pixel(0)-vc_means[0](0),2.0)+pow(vc_pixel(1)-vc_means[0](1),2.0)+pow(vc_pixel(2)-vc_means[0](2),2.0);
	dist2 = pow(vc_pixel(0)-vc_means[1](0),2.0)+pow(vc_pixel(1)-vc_means[1](1),2.0)+pow(vc_pixel(2)-vc_means[1](2),2.0);
	if (dist1<dist2)
		return 0;
	else
		return 1;
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
MCMC_LuvColor_Surface::MCMC_LuvColor_Surface()
:MCMC_RegionModel()
{
	complexity_prior = COMPLEXTY_COLOR_SURFACE;

	type = MODEL_LUVCOLOR_SURFACE;
	
}

MCMC_LuvColor_Surface::~MCMC_LuvColor_Surface()
{
	Free();
}

void MCMC_LuvColor_Surface::Free(void)
{	
}

// operators
MCMC_LuvColor_Surface& MCMC_LuvColor_Surface::operator=(const MCMC_LuvColor_Surface &x)
{
	if (this!=&x)
	{
		*((MCMC_RegionModel*) this) = *((MCMC_RegionModel*) &x);

		luv_mean  = x.luv_mean;
		luv_var	  = x.luv_var;
		if (luv_var.size()>0)
		{
			luv_var_inv = x.luv_var_inv;
			d_std = luv_var.DetNoZero();
		}
		l_surface = x.l_surface;
		u_surface = x.u_surface;
		v_surface = x.v_surface;
	}
	return *this;
}

void MCMC_LuvColor_Surface::SaveToFile(ofstream &strm)
{
	MCMC_RegionModel::SaveToFile(strm);
	luv_mean >> strm;
}

void MCMC_LuvColor_Surface::ReadFromFile(ifstream &strm)
{
	MCMC_RegionModel::ReadFromFile(strm);
	luv_mean << strm;
}

MCMC_Model* MCMC_LuvColor_Surface::CopyTo()
{
	MCMC_LuvColor_Surface *preturn;

	preturn = new MCMC_LuvColor_Surface();
	(*preturn) = (*this);
	
	return ((MCMC_Model*)preturn);
}


// compute mean and variance for the regions
bool MCMC_LuvColor_Surface::MeanFromLabel(const int nlabel)
{	
	int 		i,j;
	
	// compute the surface first
	// L channel
	l_surface.Initialize(*pmx_label,nlabel,pmx_taken,prespFilter->ras_L);
	l_surface.Estimate();
	// U channel
	u_surface.Initialize(*pmx_label,nlabel,pmx_taken,prespFilter->ras_U);
	u_surface.Estimate();
	// L channel
	v_surface.Initialize(*pmx_label,nlabel,pmx_taken,prespFilter->ras_V);
	v_surface.Estimate();

	luv_mean.SetRowVector(3);
	luv_mean.InitValue(0.0);
	nsize = 0;	

	// compute the major direction
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel)
			{				
				// get intensity
				luv_mean(0) += prespFilter->ras_L(j,i);
				luv_mean(1) += prespFilter->ras_U(j,i);
				luv_mean(2) += prespFilter->ras_V(j,i);

				nsize++;
			}
		}

		
		if (nsize>0)
			luv_mean = luv_mean / ((double)nsize);
			
		if (nsize>0)	return true;
		else			return false;
}

bool MCMC_LuvColor_Surface::MeanFromPixelSet()
{
	int 		i,j;
	bool	bgo,bgo_domain;
	PMCMC_PixelSet ppixel_set;

	// compute the surface first	
	luv_mean.SetRowVector(3);
	luv_mean.InitValue(0.0);
	nsize = 0;	

	assert(plist_domain != NULL);

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					luv_mean(0) += prespFilter->ras_L(j,i);
					luv_mean(1) += prespFilter->ras_U(j,i);
					luv_mean(2) += prespFilter->ras_V(j,i);

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

	if (nsize>0)
		luv_mean = luv_mean / ((double)nsize);

	if (nsize>0)	return true;
	else			return false;

}

bool MCMC_LuvColor_Surface::VarianceFromLabel(const int nlabel)
{
	int 	i,j;
	int			n,m;
	double		l_onsurface,u_onsurface,v_onsurface;
	double		l_pixel,u_pixel,v_pixel,value_temp[3];

	// initialization
	nsize = 0;
    
	luv_var.SetDimension(3,3);
	luv_var.InitValue(0.0);
	

	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) && (*pmx_label)(j,i)==nlabel )
			{				
				// values in l,u,v
				l_pixel = prespFilter->ras_L(j,i);
				u_pixel = prespFilter->ras_U(j,i);
				v_pixel = prespFilter->ras_V(j,i);
				// surface value
				l_onsurface = l_surface.SurfaceValue(j,i);
				u_onsurface = u_surface.SurfaceValue(j,i);
				v_onsurface = v_surface.SurfaceValue(j,i);

				value_temp[0] = l_onsurface-l_pixel;
				value_temp[1] = u_onsurface-u_pixel;
				value_temp[2] = v_onsurface-v_pixel;

				// variances
				for (n=0; n<3; n++)
					for (m=0; m<3; m++)
						luv_var(n,m) += value_temp[n]*value_temp[m];

				nsize++;
			}
		}

	
		if (nsize>1)
			luv_var = luv_var/((double)(nsize-1));

		RobustVariance(luv_var);			

		return true;
}

bool MCMC_LuvColor_Surface::VarianceFromPixelSet()
{
	int 	i,j;
	int			n,m;
	double		l_onsurface,u_onsurface,v_onsurface;
	double		l_pixel,u_pixel,v_pixel,value_temp[3];
	bool	bgo,bgo_domain;
	PMCMC_PixelSet ppixel_set;

	// initialization
	nsize = 0;

	luv_var.SetDimension(3,3);
	luv_var.InitValue(0.0);

	assert(plist_domain != NULL);

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// values in l,u,v
					l_pixel = prespFilter->ras_L(j,i);
					u_pixel = prespFilter->ras_U(j,i);
					v_pixel = prespFilter->ras_V(j,i);
					// surface value
					l_onsurface = l_surface.SurfaceValue(j,i);
					u_onsurface = u_surface.SurfaceValue(j,i);
					v_onsurface = v_surface.SurfaceValue(j,i);

					value_temp[0] = l_onsurface-l_pixel;
					value_temp[1] = u_onsurface-u_pixel;
					value_temp[2] = v_onsurface-v_pixel;

					// variances
					for (n=0; n<3; n++)
						for (m=0; m<3; m++)
							luv_var(n,m) += value_temp[n]*value_temp[m];

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

		if (nsize>1)
			luv_var = luv_var/((double)(nsize-1));

		RobustVariance(luv_var);			

		return true;
}


double MCMC_LuvColor_Surface::RobustLoglik(const double dvalue)
{
	if (dvalue>15.0)
		return 15.0;
	else
		return dvalue;
}

bool MCMC_LuvColor_Surface::LogLikelihoodFromLabel(const int nlabel)
{
	int 		i,j;
	int				m;
	double			dis,dis_uc,luvshift[3],dtmp[3];
	double			l_pixel,u_pixel,v_pixel;
	mcmcMatrix<double>	mx_variance_inv;


	nsize	= 0;
	log_pos_prob = 0.0;

	vc_vote.InitValue(0);
	crying_voice = 0.0;
	
	mx_variance_inv = luv_var.Inv();

	
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && ConsiderThisPixel(j,i) &&  (*pmx_label)(j,i)==nlabel)
			{				

				// get the values in l,u,and v
				l_pixel = prespFilter->ras_L(j,i);
				u_pixel = prespFilter->ras_U(j,i);
				v_pixel = prespFilter->ras_V(j,i);

				luvshift[0] = fabs(l_pixel-l_surface.SurfaceValue(j,i));
				luvshift[1] = fabs(u_pixel-u_surface.SurfaceValue(j,i));
				luvshift[2] = fabs(v_pixel-v_surface.SurfaceValue(j,i));

				//dis = I_u.T()*mx_variance_inv*I_u/2.0; to speed up
				for (m=0; m<3; m++)
					dtmp[m] = luvshift[0]*mx_variance_inv(0,m) + luvshift[1]*mx_variance_inv(1,m)+
					luvshift[2]*mx_variance_inv(2,m);

				dis = (dtmp[0]*luvshift[0]+dtmp[1]*luvshift[1]+dtmp[2]*luvshift[2])/2.0;

				log_pos_prob += RobustLoglik(dis);

				// compute the crying voice of this region
				dis_uc = RobustCryingVoice(RobustLoglik(dis));
				crying_voice += dis_uc;

				nsize++;
			}
		}
		
		// the third direction is treated as an independent variable
		log_pos_prob += log(sqrt(2.0*PI))*3.0*nsize;
		log_pos_prob += log(luv_var.DetNoZero())*0.5*nsize;

		// penalize regions with few valid pixels
		if (nsize <= MINIMAL_SIZE)
			log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

		log_lik = log_pos_prob;

		log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor+(double)nsize*Complexity();

		crying_voice = pow(crying_voice,0.6);			
		
		return true;
}

bool MCMC_LuvColor_Surface::LogLikelihoodFromPixelSet()
{
	int 		i,j;
	int				m;
	double			dis,dis_uc,luvshift[3],dtmp[3];
	double			l_pixel,u_pixel,v_pixel;
	mcmcMatrix<double>	mx_variance_inv;
	bool	bgo,bgo_domain;
	PMCMC_PixelSet ppixel_set;


	nsize	= 0;
	log_pos_prob = 0.0;

	vc_vote.InitValue(0);
	crying_voice = 0.0;

	mx_variance_inv = luv_var.Inv();


	assert(plist_domain != NULL);

	bgo_domain = plist_domain->GetFirst(ppixel_set);
	while (bgo_domain)
	{
		if (ppixel_set != NULL)
		{
			bgo = ppixel_set->GetFirst(j, i);
			while (bgo)
			{
				if (ConsiderThisPixel(j,i))
				{
					// get the values in l,u,and v
					l_pixel = prespFilter->ras_L(j,i);
					u_pixel = prespFilter->ras_U(j,i);
					v_pixel = prespFilter->ras_V(j,i);

					luvshift[0] = fabs(l_pixel-l_surface.SurfaceValue(j,i));
					luvshift[1] = fabs(u_pixel-u_surface.SurfaceValue(j,i));
					luvshift[2] = fabs(v_pixel-v_surface.SurfaceValue(j,i));

					//dis = I_u.T()*mx_variance_inv*I_u/2.0; to speed up
					for (m=0; m<3; m++)
						dtmp[m] = luvshift[0]*mx_variance_inv(0,m) + luvshift[1]*mx_variance_inv(1,m)+
						luvshift[2]*mx_variance_inv(2,m);

					dis = (dtmp[0]*luvshift[0]+dtmp[1]*luvshift[1]+dtmp[2]*luvshift[2])/2.0;

					log_pos_prob += RobustLoglik(dis);

					// compute the crying voice of this region
					dis_uc = RobustCryingVoice(RobustLoglik(dis));
					crying_voice += dis_uc;

					nsize++;
				}
				bgo = ppixel_set->GetNext(j, i);
			}
		}
		bgo_domain = plist_domain->GetNext(ppixel_set);
	}

		// the third direction is treated as an independent variable
		log_pos_prob += log(sqrt(2.0*PI))*3.0*nsize;
		log_pos_prob += log(luv_var.DetNoZero())*0.5*nsize;

		// penalize regions with few valid pixels
		if (nsize <= MINIMAL_SIZE)
			log_pos_prob += PENALTY_REGIONS_FEW_PIXELS;

		log_lik = log_pos_prob;

		log_pos_prob += pow((double)nsize,COMPLEXITY_POW)*scale_factor+(double)nsize*Complexity();

		crying_voice = pow(crying_voice,0.6);			

		return true;
}

void MCMC_LuvColor_Surface::GetDirectValue(mcmcVector<double> &vc, const int row_pixel,const int col_pixel)
{
	vc(0) = prespFilter->image_original.L(row_pixel,col_pixel);
	vc(1) = prespFilter->image_original.U(row_pixel,col_pixel);
	vc(2) = prespFilter->image_original.V(row_pixel,col_pixel);
}

double MCMC_LuvColor_Surface::GetDirectDis(mcmcVector<double> &vc_pixel_luv)
{
	double dtemp=0.0;

	for (int i=0; i<vc_pixel_luv.size(); i++)
		dtemp += pow(vc_pixel_luv(i)-luv_mean(i),2.0);
	return dtemp;
}


void MCMC_LuvColor_Surface::DiffuseComputedMean(const int nlabel)
{
	// compute mean
	Mean(1);
}

double MCMC_LuvColor_Surface::GetPdf(void)
{
	return 1.0;
}

// Preparation for synthesizing an image
// compute eigen vectors and eigen values
void MCMC_LuvColor_Surface::PrepreForSynthesize(void)
{
	// for the independent l channel
	mcmcMatrix<double> mx_eigenvalues;

	luv_var.Eigen(mx_eigenvalues, mx_eigehvectors);
	rf_l.Set(0.0,sqrt(fabs(mx_eigenvalues(0))),1);
	rf_u.Set(0.0,sqrt(fabs(mx_eigenvalues(1))),1);
	rf_v.Set(0.0,sqrt(fabs(mx_eigenvalues(2))),1);
}

bool MCMC_LuvColor_Surface::SynthesizeSampling(mcmcVector<double> &vc_a, const int row, const int col)
{
	mcmcVector<double> vc_temp;

	vc_temp.SetRowVector(3);

	// l,u, and v
	vc_temp(0) = (int)rf_l.Sampling();
	vc_temp(1) = (int)rf_u.Sampling();
	vc_temp(2) = (int)rf_v.Sampling();
	//vc_temp = mx_eigehvectors*vc_temp;

	mcmcMatrix<double> vc_temp1;
	vc_temp1.SetDimension(vc_temp.cols(), vc_temp.rows());
	vc_temp1 = mx_eigehvectors*(vc_temp.T()); 
	CopyMatrix(vc_temp,vc_temp1.T());


	vc_a(0) = vc_temp(0)+l_surface.SurfaceValue(row,col);
	vc_a(1) = vc_temp(1)+u_surface.SurfaceValue(row,col);
	vc_a(2) = vc_temp(2)+v_surface.SurfaceValue(row,col);

	return false; // in Luv space
}

bool MCMC_LuvColor_Surface::Reconstruct(mcmcVector<double> &vc_a, const int row, const int col)
{
	vc_a(0) = l_surface.SurfaceValue(row,col);
	vc_a(1) = u_surface.SurfaceValue(row,col);
	vc_a(2) = v_surface.SurfaceValue(row,col);

	return false;	// in Luv space
}


// get the log of the posterior probability at (row,col)
//  if pixel specified by (row,col) is already taken by other layers as marked in pmx_taken
//   then the probability should follow the law of Brownian Motion
double MCMC_LuvColor_Surface::GetLogPostProb(const int row, const int col, mcmcMatrix<int> *pmx_occupied)
{
	mcmcVector<double>	dis,luv_pixel(3);
	double			std_temp;

	assert(prespFilter!=NULL);

	if ( !prespFilter->Valid(row, col) ) return 0.0; // add by wtf

	luv_pixel(0) = prespFilter->ras_L(row,col)-l_surface.SurfaceValue(row,col);
	luv_pixel(1) = prespFilter->ras_U(row,col)-u_surface.SurfaceValue(row,col);
	luv_pixel(2) = prespFilter->ras_V(row,col)-v_surface.SurfaceValue(row,col);

	// u and v
	if (d_std >= 50.0)
	{
		dis = luv_pixel*luv_var_inv*(luv_pixel.T())/2.0;
		std_temp = d_std;
	}
	else
	{
		// if variance is too small then the force will be falsified
		mcmcMatrix<double> mx_temp;
		mx_temp = luv_var_inv/pow(100.0/d_std,1.0/((double) luv_var_inv.rows()));
		dis = luv_pixel*mx_temp*(luv_pixel.T())/2.0;
		std_temp = 50.0;
	}

	return	log(sqrt(2.0*PI))*3.0+log(std_temp)*0.5*2.0+dis.Data(0);
}


void MCMC_LuvColor_Surface::OutputInfo(char *info)
{
	if (info!=NULL)
		sprintf(info, "%.2f(%.5f) surface %.1f:%.1f:%.1f--%.1f", log_pos_prob,log_lik/nsize,
		luv_mean(0),luv_mean(1),luv_mean(2),d_std);
}
